next up previous contents
Next: 4. CASA Imager Memory Up: Imaging Algorithms in CASA Previous: 2. Currently available Imaging   Contents

Subsections


3. Imager Refactor

3.1 Requirements

Provide the following use cases from a consistent code base :

  1. Data Selection : Selection from a single MS, a list of MSs, or a MMS. Ability to select 'data' vs 'corrected_data' for imaging. Add weight-scaling options. Add readonly mode, in addition to a choice between virtual or scratch model columns.

  2. Image definition : Multi-frequency-synthesis (nterms=1 or more) and Cube (parameters specifiable in multiple spectral frames and index/frequency/velocity conventions).

  3. Deconvolution options : Hogbom, Clark, ClarkStokes, Multi-Scale, Multi-Term MFS. Re-use existing numerical code.

  4. Imaging/Gridding options : Standard gridding (prolate spheroidal), W-Projection, A-Projection with convolution functions from Airy disk beams, ray-traced beams externally supplied beams, and with and without squint correction. Re-use existing numerical code where possible. All modes must work for mosaics (joint residual images) and multi-term MFS.

  5. Outlier image fields : Combinations of all of the above must work for one or more image fields. Iteration and mask control for multi-field images must be consistent across algorithms. Allow different imaging and deconvolution parameters for different fields.

  6. Mask generation : For all above options, masks to be constructed from region specifications, external images (in any coordinate system) with either T/F or 1/0 masks, interactively drawn masks, autoboxing. Allow mask merging (e.g. autoboxing within a given pb level) and mask extensions (across channels).

  7. Iteration control : All algorithms must respond in the same way to iteration parameters ( thresholds, scale-factors, gain, and iteration counters across image fields, channels and stokes planes, major-cycle triggers, mask-drawing triggers etc... ). Perhaps allow field and channel-dependent controls.

  8. Runtime interactivity : Allow some iteration control parameters to be modified at runtime. In particular, a 'stop' option to terminate minor cycles.

  9. Parallelization : Must allow all of the above features and interaction with parallelized major cycles for MFS, and full parallelization (major and minor) for Cube imaging. All parallelization must be implemented via Python. Iteration control must behave the same for serial and parallel runs (where logically possible). Serial and parallel runs must use the same underlying numerical code and produce the same results.

  10. Log messages : Results of setup options (e.g. number of unflagged visibilities, psf beam size, theoretical image RMS, etc ) must be shown in the log. The progress of deconvolution (e.g. peak residual, iteration count, relation to thresholds) must be shown in easily readable format ( especially for multi-field, many-channel cube runs ).

  11. Multiple tasks : One large 'clean' task that does everything. Perhaps also different flavours of the 'clean' task to expose parameters by complexity level. All these tasks are meant to replace the current 'clean','ft','setjy','deconvolve','pclean' tasks.

  12. Smaller tasks : Provide smaller tasks to run each step of imaging and deconvolution separately ( 'makepsf', 'makeresidual', 'deconvolve','setmodel' ) ensuring that they construct and write out only required images. All these tasks should use the same underlying code.

3.2 Software Design

These slides from 2013 CASA Developer Meeting contain an outline of the required use cases, the modules being built in the refactored Imager classes and how they connect together for serial and parallel runs.

3.2.1 C++ Classes

  1. Image Sets : ImageStore is a C++ class that holds the residual, model, psf, restored, mask, weight(pb) images for one field. This is the image unit that is constructed and passed around between other classes that operate on it. For now, all transfers of images between modules (and across compute nodes for a parallel run) is done via disk, with each module constructing its own reference to the set of images on disk.

  2. Major cycle : SynthesisImager is the C++ interface layer that controls data selection from one or more measurement sets, and the definition and creation of images. It maintains one or more mappings of gridding methods to images (single/multi-field). Gridding methods are implemented with FTMachines. SIMapper is a wrapper that maps one or more FTMachines to one SIImageStore, and SIMapperCollection is a list of such mappers used for multiple fields. The SynthesisImager contains one or more MeasurementSets, MSSelections, and one VisibilityIterator for data access, and passes only VisBuffers to the Mappers.

  3. Minor Cycle : SynthesisDeconvolver is the C++ interface layer that controls deconvolution algorithm definition and execution, for one field. The SDAlgorithmBase class implements loops through channels and stokes planes, calls virtual methods for actual iteration steps, and checks a common iteration controller (SIMinorCycleController) for stopping criteria. Derived SDAlgorithmXXX classes implement the numerical deconvolution algorithms, and through their base class, get identical iteration and channel/stokes plane controls.

  4. Iteration control : SIIterBot is the C++ interface layer to iteration and interaction control. Thresholds, iteration counts and major-cycle triggers are implemented here in one place, for all minor cycle algorithms to use. The SynthesisDeconvolver takes in a subset of the SIIterbot, the SIMinorCycleController, and reads and updates its iteration control and imaging summary parameters. The SIIterBot is synchronized at major/minor cycle boundaries and all deconvolvers see synchronized parameters.

  5. Normalizations : SynthesisParSync is the C++ interface layer to a class that handles pre-deconvolution and pre-degridding normalizations of images by weight images. In the simplest case, this is simply a division by sum of weights before deconvolution. For wide-field imaging, it controls flat-sky vs flat-noise normalizations in both directions. For parallel runs, it allows for a gather step before normalization.

  6. Tool interfaces : One-to-one mappings were provided between the C++ interface layers and python, with no extra functionality implemented at the cmpt level. The tools created are synthesisimager, synthesisdeconvolver, synthesisparsync, synthesisimstore and synthesisiterbot.

  7. Error Handling and Exception messages : Exception try/catch clauses at all major levels, with messages pre-pended to the exception message at each level to help identify sources of exceptions in the future.

Descriptions to be added when implemented : Mask handling, Runtime Interaction.

3.2.2 Python Classes

The SynthesisImager contains a SIMapperCollection to map one set of data to multiple image fields. The SynthesisDeconvolver operates on one field. The connection between the two is maintained one level above these classes, currently in Python through the individual tool interfaces. Note however, that all python functionality can be duplicated within C++ if later desired. The current choice of Python was driven by the requirement that all parallelization be done via python.

  1. Wrapper Class : PySynthesisImager is a base class that strings together all relevant tools to implement broad functions such as initializeImagers/Deconvolvers/ParallelSync/terationControl, makePSF, runMajorCycle, runMinorCycle, runMajorMinorLoops, getSummary, and deleteTools. This class can hold one synthesisimager, and one or more synthesisdeconvolvers (for multi-field imaging) and synchronizes their iteration parameters via a common synthesisiterbot. By choosing the appropriate function calls, major cycles, psf generation, and minor cycles can be triggered in isolation.

  2. Parameter Handling : ImagerParameters is a class that does parameter verification and type-checking for each functional block (data selection, image definition, gridding option, deconvolution options, iteration options) before passing it into C++. All allowed and forbidden combinations (as functionality is internally added) are controlled here. Several old and new field-specific imaging and deconvolution parameters can be specified in an input outlier file.

3.2.3 Parallelization

All parallelization is currently done at the python level. The basic wrapper modules have been adapted to run parallel major or minor cycles or both.

  1. Cluster management : PyParallelImagerHelper is a class that manages a compute cluster (via simple_cluster) and provides interfaces to start it up, send commands to casapy shells on remote nodes, check for job completion, and shut it down. It also implements very basic heuristics for data and image partitioning.

  2. Parallel major cycles : PyParallelConSynthesisImager is a derivative of PySynthesisImager which parallelizes the major cycle. Instead of holding a single synthesisimager tool, it creates multiple synthesisimager tools, one on each available node. Calls to the tool methods are sent via PyParallelImagerHelper.

  3. Parallel minor cycles : PyParallelDeconvolver is a derivative of PySynthesisImager which parallelizes the minor cycle on image fields. Instead of holding one or more synthesisdeconvolvers, it creates one tool on each available node. This is a simplified form of parallelization that may not be very useful.

  4. Cube parallelization : PyParallelCubeSynthesisImager is a class that uses PyParallelImagerHelper to deploy multiple instances of PySynthesisImager, one on each available node. This does independent imaging and deconvolution on each frequency chunk, with no interaction between the chunks for iteration control. There is currently no concatenation of the image cube at the end.

    An alternate approach to parallel cube imaging can be implemented as a combination of PyParallelContSynthesisImager and PyParallelDeconvolver, where the major and minor cycles are separately parallelized and iteration control is synchronized across the frequency chunks via a common synthesisiterbot.

3.3 Current status

3.3.1 Work completed as of October 2013

C++ and Python classes have been written as per the above code design, to support the following :

Multiple image fields Working, with field-dependent parameters specified via an outlier file
Cube Imaging Working, for spectral coordinates specified in frequency
MFS Imaging ( nterms=1 ) Working
Multi-term MFS Imaging ( nterms>1 ) No
Deconvolution options Hogbom Clean
Gridding options Standard (prolate spheroidal) and W-Projection
Mosaic Possible with standard gridding. No primary beams yet.
Centralized Iteration Control Working for niter, cycleniter, threshold, cyclethreshold, gain
Centralized minor cycle loops over field, channel, stokes planes Working
Major cycle parallelization (for MFS and Cube) Working
Cube Parallelization Working, without iteration synchronization between major and minor cycles. A version with synchronization is possible but not implemented
Minor cycle parallelization (for multi-field) Working, but not intended for use
Mask Handling (input masks,interactive masks, autoboxing) No
Runtime parameter modification Prototyping in progress, via a new pop-out Viewer GUI.
Cube specification in all spectral frames No
Clear log messages Trying. Experimenting with a summary plot to track deconvolution progress
Input parameter validation Basics are working in python wrapper layer. Needs more
Tool Interfaces Done
Python Wrappers Basics done for single/multi-field, serial/parallel. Needs more details
Task Interface No
Test programs Python function level tests for all supported cases

3.3.2 Remaining work

  1. Test/finish current functionality

    UR : Convert all existing tests into proper functional tests for all available use cases (single/multi-field, cube/mfs, serial/parallel major cycles, full-clean/major-only/minor-only). Add python level tests for faceted imaging.

    UR : Add the image restoration step to SDAlgorithmBase, and provide a python level function for this.

    UR : Provide Python wrapper functions for partial operations : makePSF, makeResidual, predictModel, etc.

    UR/KG : Memory usage tests to evaluate the current approach of each module constructing SIImageStore objects from disk. If it is not good enough, then either focus on implementing live transport of SIImageStores for serial runs, using the prototyped transfer of tool instances across the C++/Python boundaries, or postpone this optimization to when parallelization gets away from Python.

    KG : Fix some data selection bugs in SynthesisImager ( multiple SPWs, multiple MSs )

    KG : Write data and image partitioning code in PyParallelImagerHelper.

  2. Task interface

    UR : Write task interface and connect to top level Python classes.

  3. Spectral Frames

    TT : Implement a function in synthesisimager_cmpt.cc or SynthesisImager to cover all the spectral frame specification use cases as documented.

    TT : Write functional tests for all spectral-frame use cases ( single field, cube ), making only dirty images.

  4. Mask Management

    UR : Provide an interface to connect SynthesisDeconvolver to a basic MaskHandler

    TT : Implement insides of MaskHandler to make masks from input region specifications, input images, PB levels, autoboxing, and merging of multiple types of masks. The interface would be parameters, plus image names ( or imageinterfaces ) for the residual image, psf, mask, pb. Minor cycle algorithms would use this class, or derived versions of it ( for example, multi-scale expands masks, multi-term spreads masks across Taylor terms ).

    TT : Evaluate the feasibility of replacing the insides of the current 'makemask' task with this MaskHandler ( via a new maskhandler tool interface ).

  5. Runtime Interactivity

    DS : Continue on interactive GUI prototype, to edit parameters like 'niter', 'threshold' or 'gain' and have this modify parameters directly in the SIIterBot class, from which all SIMinorCycleControllers will get updated parameters at their next checkpoint.

    UR : Test that the iteration controller responds appropriately to these runtime changes. Check performance trade-offs between full interactivity (check for changes after each iterations) vs only at predefined checkpoints (set via the SIIterBot).

    UR : Evaluate runtime parameter modification behaviour with parallel minor cycle runs, and decide what to do.

    DS/UR : Connect interactive GUI to MaskHandler, sending region info and autoboxing parameters into the MaskHandler and asking for, receiving, and rendering an updated mask. This communication would be through the SIIterBot, which would trigger a mask-update step from SDAlgorithmBase.

  6. Setjy Functionality

    UR/KG : Provide a path for 'setjy' functionality to be added into the new framework and run it from python.

    TT/KG : Connect/test setjy functionality, and make functional tests.

    TT : Make new setjy task interface using PySynthesisImager (after the clean task interface is done).

  7. Other deconvolution algorithms and gridders

    UR : Connect other standard minor cycles to new framework ( clark, clarkstokes, multiscale )

    UR/KG : Connect multi-term wrapper FT-machine, and MS-MFS minor cycle algorithm.

    KG/SB : Connect A-Projection gridders to get proper mosaics.

3.4 Relevant design/discussion documents


3.4.1 New Clean task parameters

This section lists an alternate tasking interface for image-reconstruction in CASA. This interface follows the underlying structure of the algorithms being used, and splits the parameters into five sections : data selection, image-definition, iteration-control, imaging, deconvolution.

It will be possible to run only 'selection + image-definition', generate some logger output ( Number of visibility points for image-rms calculations, selected frequency-ranges in the data-frame and LSRK, etc.. ) and then exit. This is to allow verification of image-coordinates before proceeding with the full imaging.

A run with niter=0 would do one major cycle to make PSFs and residual images (and return a dictionary of info). If a modelimage is specified, this step will save or predict the model (the equivalent of task 'ft' but with full flexibility of setup and ftmachine choices). Irrespective of 'niter', this first major cycle will generate logger information about sum-of-weights, expected image-rms, min and max spatial-scales spanned by the uv-coverage, clean-beam size, PSF sidelobe level, etc. If possible, if no parameters are changed, it should be possible to re-start deconvolution directly from the minor cycles, and avoid re-doing the major cycle.

3.4.1.1 Data Selection

Select a subset of a dataset (or multiple datasets). All selected data will be used for image reconstruction.

-> Imager tool interface : im.open(); im.selectvis();

vis '' or ['',''] Name of input visibility file (for multiple ms, give a list of strings)
field '' or ['',''] Field Name or id (for multiple ms, give a list of strings)
spw '' or ['',''] Spectral windows e.g. '0$ \sim$ 3', '' is all (for multiple ms, give a list of strings)
selectdata True More selection (for multiple ms, give a list of strings)
timerange '' or ['',''] Range of time to select from data
uvrange '' or ['',''] Select data within uvrange
antenna '' or ['',''] Select data based on antenna/baseline
scan '' or ['',''] Scan number range
  False No more data selection

Keep the 'selectdata' parameter even though its default is True and it doesn't help save screen-real-estate. Jeff has a plan to consolidate 'selection' across tasks, so until then, don't change anything.

3.4.1.2 Image Definition

Define the coordinate systems of output images.

-> Imager tool interface : im.defineimage();

Postpone discussion about mode='cube' vs mode='channel','velocity','frequency', as the CSSC does not yet agree on the 'cube' idea/document circulated by JO. For now, no changes here.

The following is a list of points to discuss, based on suggestions from KG and DP.

  1. Remove 'outframe' because Imager always grids in LSRK, and provide an external tool/task to 'set' the conversion-layer in the image-coordinate system for use by the viewer and fits-writers for on-the-fly conversion.
  2. Alternatively, add support for 'native-gridding' vs 'LSRK' gridding, to support datasets that are pre-gridded to the target image-channelization, and which do not need/want the LSRK conversion.
  3. 'reffreq' is decided purely by data-selection (middle of the range in LSRK)
  4. restfreq' : Needed before conversion to FITS. What about multiple rest-freqs ?
  5. Perhaps remove support for velocity-image-coords, and 'veltype', and force this step as a post-deconvolution conversion ?
  6. Add a new parameter 'excludefreqs' either here, or in data-selection. It would be a list of frequency ranges (or velocities), with units and ref-frames, to exclude from gridding/imaging. This is to prevent the user having to guess/calculate which TOPO data channels would map to these LSRK frequencies or velocity-ranges at at-least one timestep.

Postpone decisions about the existence of parameters 'reffreq', 'restfreq', 'outframe', until after KG's suggestions for reducing user-confusion about spw-selections vs image-freq-definitions, especially when frequency-frame changes require frequency-regridding, and there is NO one-to-one mapping between data channels and image channels.

mode 'mfs' Combine all frequencies during imaging to produce a wideband image. The bandwidth is computed from the selected data-range, converted to LSRK.
reffreq ' ' Reference frequency of the continuum image. By default, it is the middle of the frequency-range computed from the data-range converted to LSRK. This is also the ref-freq used by the 'msmfs' deconvolution algorithm for its Taylor-polynomial model of the source spectrum.
nterms 2 Number of Taylor-polynomial terms used to model the sky spectrum.
  'channel' Cube, image-channels are defined via data-channel-frequencies converted to LSRK at first timestep.
nchan -1 nchan
start '' start chan (from zero in the MS, ignoring data-selection)
width 1 step size
interpolation 'linear' Interpolation to use when regridding visibilities from data to image-frequencies. As a default, this should be 'nearest' for mode='channel', but we currently cannot have different subparam defaults for different parent-params
outframe 'LSRK' The image is always made in LSRK. An external frame-convertor can be set to one of the following 'LSRK' 'LSRD' 'BARY' 'GEO' 'TOPO' 'GALACTO' 'LGROUP'.
  'frequency' Cube, with image-channels defined via frequency-ranges.
nchan -1 nchan
start '' start freq
width 1 step size
interpolation 'linear' Interpolation to use when regridding visibilities from data to image-frequencies.
outframe 'LSRK' The image is always made in LSRK. An external frame-convertor can be set to one of the following 'LSRK' 'LSRD' 'BARY' 'GEO' 'TOPO' 'GALACTO' 'LGROUP'.
  'velocity' Cube, with image-channels defined via velocity-ranges.
nchan -1 nchan
start '' start velocity
width 1 step size
interpolation 'linear' Interpolation to use when regridding visibilities from data to image-frequencies.
outframe 'LSRK' The image is always made in LSRK. An external frame-convertor can be set to one of the following 'LSRK' 'LSRD' 'BARY' 'GEO' 'TOPO' 'GALACTO' 'LGROUP'.
veltype 'radio' Velocity definition of output image. 'optical' or 'radio'.
stokes 'I' Stokes params or correlation products to image : I,Q,U,V, RR, LL, XX,YY and all combinations allowed by the input data (I,Q,U,V,IQ,IV,QU,UV,IQU,IUV,IQUV,RR,LL,RRLL,XX,YY,XXYY).
imagename ' ' Pre-name of output images ( [ ' ', ' ' ] for multi-fields )
imsize [256, 256] x and y image size in pixels. Single value: same for both. ( [ [,], [,] ] for multi-fields )
cellsize ['1.0arcsec', '1.0arcsec'] x and y cell size(s). Default unit arcsec. ( [ [' ',' '],[' ',','] ] for multi-fields)
phasecenter ' ' Image center: direction or field index ( [ ' ', ' '] for multi-fields)
mask ' ' or [ , , ] Cleanbox(es), mask image(s), region(s), or a PB-level ( [ [ , , ] , [ , , ] ] for multi-fields)
modelimage ' ' or [' ' , ' '] Name of model image(s) to initialize cleaning ( [ [ , , ] , [ , , ] ] for multi-fields)
outlierfile ' ' Outlier file (new format) with imagename, imsize, cell, phasecenter, mask, modelimage. A single main field is specified in the task-parameters, and only outliers are in this file. Outlier files cannot be mixed with list-modes.

(1) A convertor can be provided for the AIPS format, but without 'modelimage' support. i.e. no mixing of inline-lists and outlier-files for parameters un-supported by the AIPS format. If the convertor is run separately, the user can then edit the files to add model-image info.

(2) Old CASA-Imager outlier-file format will go away.

Example outlierfile format : Two outliers, supplying a mask and modelimage only for the first. These two fields are in 'addition' to the main single field specified in the task parameters. 

    imagename = 'M1_1' 
    imsize = [128,128]
    phasecenter = 'J2000 13h30m52.159 43d23m08.02'
    mask = ['out1.mask', 'circle[[40pix,40pix],5pix]' ]
    modelimage = 'out1.model'

    imagename = 'M1_2' 
    imsize = [128,128]
    phasecenter = 'J2000 13h24m08.16 43d09m48.0'

Inline lists for multi-field specification : Allow only for simple cases, but force outlier files as the only option for modes that require multiple levels of nested lists, which can also be different levels for different parameters (for example, multi-term + multi-field + multiple-masks per field).

Initial implementation will allow only outlier files and no inline-lists. Support for inline lists will then be added for 'specific, simple' cases for ease-of-use.

3.4.1.3 Iteration Control

Choose the type of iterations (major and minor cycles), gain, and thresholds. All iteration-counts are maintained per-field, per-channel and per-stokes.

-> im.clean()

niter 160 Maximum number of minor-cycle iterations, counted across major-cycles. This is over-ridden by the user-specified stopping-threshold.
npercycle -1 Force hogbom-clean, with major-cycles only at the beginning and end. Do not use FluxLimit.
  50 Maximum number of minor-cycle iterations before triggering a major-cycle. This is overridden by 'FluxLimit' calculations controlled by sub-parameters, to trigger a major cycle. Also overridden by the user-specified stopping-threshold.

Total number of major cycles $ >=$ niter/npercycle

For example, niter=160, npercycle=50. Due to the flux-limit, an extra major-cycle is triggered after iteration 28. Major-cycles will occur after iterations 0, 28, 78, 128, 160.

FluxLimit = peakResidual x fractionOfPsf where fractionOfPsf = MIN(maxPeakFraction, MAX( minPeakFraction, PSFsidelobe x cyclefactor) )

This FluxLimit is controlled by the PSF sidelobe level and 'cyclefactor'

maxpeakfraction 0.8 The highest fraction of the peak residual that FluxLimit can be. Must be smaller than 1.0 !
minpeakfraction 0.1 The lowest fraction of the peak residual that FluxLimit can be. Must be larger than 0.0
cyclefactor 1.0 Scale-factor for the psf-sidelobe level. Rename to 'sidelobescalefactor' ?
nitertol 0.05 If we get to npercycle iterations, or a fluxlimit, and less than 5% of total iterations remain, stop and tell the user those last-few iterations won't be done. If set to 'zero', then those last few iterations will be done after the major cycle. May be dangerous to just do those extra iterations as it may go too deep.
checkdivergence 0.5 If the peak residual increases by this fraction, do not use the last model, but trigger a major cycle
negcomponent -1 If larger than zero, stop minor cycle iterations when this many negative components are found. When 'multiscale' is turned on, this count applies to negative components at the largest spatial scale.

gain 0.1 Loop gain for cleaning
  -1 Adaptive loop gain, following the Levenger-Macquart method, to go faster when possible, and slow down when in a shallow part of chi-square space (not discussed yet, but already available/tested for msmfs in current task)
threshold '0.1mJy' Main stopping threshold for each field, channel and stokes-plane. If any field, channel, or stokes-plane reaches this threshold, stop minor cycles for the current field, channel or stokes-plane, and move on to the next field, channel or stokes-plane. This overrides all other iteration-control.
interactive True Pop up a GUI for mask-drawing with (stop, continue-till-end, continue until next stopping-point). Stop after every major cycle. Also allow modification of niter, npercycle, maxpeakfraction, minpeakfraction, cyclefactor, threshold, gain. This will be implemented via a text-file on disk, but will have a GUI (for ease of use).

If ftmachine='refpb' or 'aproject', always overlay a contour of the PB at 'pblimit'. If pbcor = True, overlay a contour of the PB at 'pbmin', in a different colour.

If calready=True, save the model at every stopping-point, so that one can exit at any time.

pauseoption False Stop only at major-cycle boundaries
  True A button on the interactive-clean GUI stays alive, and can be clicked to 'pause', even between minor-cycle iterations. Then, masks can be modified and all those parameters can also be modified via the GUI. An exit from with minor-cycle iterations, will trigger one last major cycle.
plotprogress False Nothing
  True Show continuously-updating plots of total flux and peak residual

  False Run without a GUI, but allow runtime parameter modification via a text file on disk containing parameters. This will be the same text-file involved in parameter-modification via the GUI. File name will be [imagename].control. Decide later if 'plotprogress' should be possible here too.

To start with, do not implement `` chaniter : Clean each channel (minor and major cycles) completely before proceeding to the next ``. If later required (purely as a usage-mode and not for performance), it can be implemented at the top level, with the expected caveat of diminished performance, because of (a) ignoring data-locality in tiles, and (b) having to read multiple channels even to make one image channel because of freq-frame conversions.

An option to do a flat-sky normalization before showing the current residual image was requested. It is possible to do this extra normalization, but to avoid making an extra copy of the whole image, it would mean an extra sequence of divisions/multiplication. This is an extra step, and with small numbers near the edge, may lead to numerical instability. Therefore, the initial implementation will only 'show' PB-contours at the 'pblimit' level when ftmachine='refpb' and 'aproject'. The extra division/multiplication by the PB can be added later if still required.

3.4.1.4 Imaging Options

Choose weighting schemes and direction-dependent effects to apply while gridding visibilities and constructing the image to be sent to the deconvolution algorithm.

-> Imager tool interface : im.weight();im.setvp();im.setoptions();

Note : Mosaicing (phase-gradients on gridding-convolution functions) is implicit in all of the following, and is automatically triggered when multiple pointings are present in the selected data (i.e. no parameters required). To do a linear-mosaic with separate deconvolutions per facet, use multi-field options followed by im.linearmosaic()).

weighting 'natural' Natural weighting
  'uniform' Uniform weighting
  'briggs' Robust weighting
robust 0.0 go between uniform and natural
npixels 0 npixels to decide uv-cell size over which weights are computed
  'briggsabs' Robust weighting - 2
robust 0.0 go between uniform and natural
noise '1.0Jy' noise estimate
npixels 0 npixels to decide uv-cell size over which weights are computed
uvtaper False No uv-tapering
  True Apply additional uv tapering of visibilities
outertaper [] uv-taper on outer baselines in uv-plane
innertaper [] uv-taper in center of uv-plane ?? (not implemented)
ftmachine 'ft' Use the standard prolate-spheroidal for gridding.
  'sd' Use the FT of the single-dish PB as the gridding-convolution kernel
  'refpb' Use the FT of antenna-dependent PBs, at ref-freq, and one time as the gridding-convolution kernel current ftmachine=mosaic. Data weights are not considered in PB computation. Normalization for the minor-cycle is always flat-noise.
pblimit 0.02 limit beyond which to ignore
  'aproject' Use time-varying, frequency-dependent, and antenna-dependent gridding convolution functions constructed as FTs of PBs, rotated and scaled for time and freq variations, squinted, and weighted according to data-weights. Normalization for the minor cycle is always flat-noise..
pblimit 0.02 limit beyond which to ignore
cfcache 'sss' cache directory

wprojplanes

 1 No w-projection. Use the gridding convolution function selected via 'ftmachine'.
  64 Use w-projection kernel (FT of Fresnel propagator) with this-many w-planes, in addition to whatever selected by 'ftmachine'

facets

 1 For each image-field, use a set of 'facets x facets' sub-images during gridding. Deconvolution will always see a joint image (covering the area : 'facets x facets'). Faceting can be used as an alternative to w-projection (wprojplanes=1), along-with W-projection (wprojplanes$ >$ 1), or along-with A-projection or 'ft'. In all cases (especially mosaicing), the convolution-functions used will calculated for the facet-size field-of-view.
pbcor False Output image represents $ I^{sky} \times PB$ . No extra normalization.
  True Output image represents $ I^{sky}$ . Do a post deconvolution correction using an average PB. If nterms$ >$ 1, this does a post-deconvolution wide-band pb-correction.
pbmin 0.1 Mask for post-deconvolution pb-correction

Initially, not all possible valid combinations will be enabled. Un-supported combinations will be checked-for and stopped at the start of the task, and enabled when the internals are ready for it. The main point here is that the interface will support all physically-valid combinations in a way that follows the mathematics of the imaging process.

3.4.1.5 Deconvolution Options

Choose the deconvolution algorithm based on sky parameterization. This will run separately per field, channel, and stokes-planes (except when combinestokes=True). Iterations are the number of minor-cycle iterations decided by iteration-control. This part will also be available at the tool-level as stand-alone image-based deconvolvers.

-> Imager tool interface : im.setscales(), im.settaylorterms(), im.clean();

deconvolver 'hogbom' A variant of the Hogbom clean algorithm. Use a point-source flux model. Use full PSF.
  'clark' A variant of the Clark-clean algorithm. Use a point-source flux model. Use truncated PSF, minor cycle includes periodic subtraction using gridded visibilities.
  'clarkstokes' Similar to 'clark', but search for peaks in $ I^2 + Q^2 + U^2 + V^2$ . Requires a 4-stokes-plan input image.
  'multiscale' Model the sky as a collection of Gaussian-like components of different sizes.
scalevector [0,3,10] List of scale sizes in units of pixels, or as strings : ['2.0arcsec', '6.0arcsec','20.0arcsec'].
smallscalebias 0.6 Give more weight to small scales
  'mem' maximum entropy (SUBPARAMS)
  'asp' asp-clean (SUBPARAMS)
autobox False No autoboxing
  True Do autoboxing. Similar to interactive=True one could contol whether the autoboxing should be repeated every N minor-cycle iterations or at major-cycle boundaries. (Add autoboxing SUBPARAMS)
restoringbeam '' or ['','',''] Gaussian restoring beam for CLEAN image. Default is to use a frequency-dependent beam, fitted from the PSF main lobe at each frequency.
calready False Do not save a model anywhere inside the MS.
  True Save a model for later use in self-calibration or plotting.
usescratch False Save model as a virtual model-column record. The ftmachine state, with the image and data-selection info are stored.
  True Write to the model-data scratch column in the MS.
async False If true the taskname must be started using clean(...)

3.4.2 Spectral Frames in the Refactored Imager/Clean

Spectral coordinate systems produced by the CASA Imager are always tied to a specific reference spectral frame (we call this the base frame). Channel labels are in units of frequency, referenced to this base frame. Channels can be either regularly spaced, or irregular. An optional conversion layer can be added on top of the base layer, which allows on-the-fly relabeling of spectral channels into different spectral frames. Access functions also exist to read out the spectral channel labels in velocity instead of frequency, also on-the-fly.

The following sections describe the spectral coordinate systems in images produced by clean in terms of the base frequency reference frame and conversion layer. For use-cases where on-the-fly relabeling via conversion layers will not suffice, further operations (outside of Imager) to explicitly change the base reference frame and possibly regrid the channels, will be provided (imregrid, imreframe, exportfits).

Section 3.4.2 describes the mathematical operations involved in mapping data channels to image channels. Section 3.4.2 describes three supported software Doppler tracking options with a list of parameters relevant to each and the spectral-frame structure of the output image. Section 3.4.2 describes how the spectral axis can be defined by the user, and how various parameters are interpreted. Section 3.4.2 discusses ways in which the output images and their conversion layers can be used for further analysis. Section 3.4.2 contains some details about frame definitions and conversions supported by the casacore package, which Imager uses.


3.4.2.1 Data channels to Image Channels

During the gridding process, a choice is made on how to map data channels to image channels, based on the channelization encoded in the MS and the user-specified image channelization. There are three reasons to deviate from a one-to-one mapping of data channels to image channels. This section applies to both the old and new Imager.

  1. If the user-specified 'start' frequency is not at a data channel boundary, the visibility data and weights are interpolated to evaluate them at the centers of the shifted frequency grid.

  2. When image channels are wider than data channels, visibilities and weights are binned and gridded via MFS within each image channel.

  3. On-the-fly software Doppler tracking can also be done as a time-dependent frequency relabeling of the data channels, before interpolation and binning choices are made. There are three options, described in section 3.4.2.

Usually, a combination of these three operations are performed during gridding. There are off-course special cases where only some or none of them apply.


3.4.2.2 Software Doppler tracking options

For the purpose of this document, a time independent frame is a frame in which the source observed has not changed its velocity over the period of time of the observation. A time dependent frame is one in which the observed source has changed its velocity during the observation period. If datasets from different epochs are combined during imaging, the relevant period of time to consider is the entire range spanned by all datasets.

There are three software Doppler tracking options, which will be controlled at the task level. Individual parameter interpretations are described in section 3.4.2.

  1. mode='cube': Convert data frequencies to the LSRK time-independent spectral frame.
    Output image base frame : LSRK
    Conversion layer : Yes

    Data frequencies are read as vb.lsrfrequency(), and compared with image channel bins also in lsr-frequency, to decide which data channels go into which image channel bins. This is a time-dependent relabeling of data frequencies and spectral lines from astrophysical sources are lined-up in this process.

    Relevant User Parameters : start, step, nchan, frame, veltype, restfreq. Internally, this is converted to LSRK frequencies to construct an output CASA image with a base frame of LSRK.

    At the end of clean, a conversion layer is set to allow easy access to the list of frequencies in the 'frame' in which the user specified the input : LSRK, LSRD, BARY, GALACTO, LGROUP, CMB. Note that for TOPO and GEO, which are actually time-dependent frames, the conversion layer encodes a specific time at which the conversion is valid. The convention in 'clean' will be the start epoch of the dataset, but this can be changed via the imreframe task with outframe='topo'.

  2. mode='cubesrc' : Convert data frequencies to the SOURCE frame.
    Output image base frame : SOURCE
    Conversion layer : No

    If the FIELD table of the source being imaged contains ephemeris information, a time-dependent relabeling of the data frequencies (software Doppler tracking) is done to make spectral lines stationary in the moving source frame.

    If the FIELD table of the source being imaged does not contain ephemeris information (i.e. the source is not a solar system object), the software Doppler tracking will follow a conversion to LSRK. In addition, a systemic velocity has to be specified w.r.to a spectral frame, which will be recorded in the image for future velocity readouts.

    Relevant User Parameters : start, step, nchan, frame, veltype, restfreq, sysvel, sysvelframe. The base frame of the output image will always be SOURCE. sysvel, sysvelframe represent the systemic velocity w.r.to a specific frame that will be embedded in the coordinate system (for future velocity readouts). These two parameters are ignored if the ephemeris information is available. This is the only mode that allows start,step to be specified in frame='SOURCE' in addition to other standard frames.

    ( Note : It is currently not planned, but if needed, another parameter can be added in this mode, to explicitly point to an external ephemeris table [ref. GM's example about observing Titan but getting RFI from Cassini !]. )

    No conversion layer will be set at the end of 'clean'.

  3. mode='cubedata' : A one to one mapping between data channels and image channels.
    Output image base frame : REST, NOFRAME.
    Conversion layer : No

    No time-dependent frequency frame transformations are done during gridding/binning. In this case data frequencies are read as vb.frequency() and compared directly with image frequencies. With a dataset observed in a time-dependent frame, astrophysical spectral lines remain smeared out. Only local signals such as RFI will stay un-smeared (terrestrial RFI in case of TOPO data)

    Relevant User Parameters : start, step, nchan, veltype, restfreq. Here frame is not an option as start, veltype, restfreq will be interpreted literally to construct the frequency grid, with no further conversions.

    If the MS frame is labeled as 'REST' (in all input MSs), the output image will also be 'REST'. mode='cubedata' is the only mode in which imaging of 'REST' MSs will be allowed.

    For all other input MS frames, the output image will be labeled with a (new) frame name of 'NOFRAME', to represent an arbitrary frame. This frame has to be completely generic ( it has to be valid for joint imaging of multiple MS's that could be in different reference frames ) and there is no logical true frame with which to label the output image.

    We will provide a helper function to change 'NOFRAME' to something else, after imaging. This function will operate only on images containing 'NOFRAME' as the base frame.

    Also note that even when there is no time-dependent relabeling of frequencies, interpolation and binning can still happen if the desired image channel boundaries do not match those in the input dataset.

    No conversion layer will be set at the end of 'clean'.

  4. mode='mfs' : Multi-frequency synthesis, where there is only one large image channel. This will always be in LSRK, with the image frequency coordinates decided by the spw data selection parameter. There will be no conversion layer.

3.4.2.3 Imaging a pre-Doppler-tracked dataset

An MS coming out of cvel or mstransform will be treated the same way as any other MS observed directly in the frame that the MS is labeled with.

  1. An mstransformed dataset labeled in a time-independent frame ( LSRK, LSRK, BARY, etc.... ) can use mode='cube'. The base frame of the output image will be in LSRK, with a conversion layer based on how the user chose to specify input parameters. If the mstransformed MS is already in a time-independent frame, the code will detect that no additional time-dependent frequency shifts are required. A similar situation holds for datasets labeled in the SOURCE frame and when mode='cubesrc' is used.

  2. Datasets that need channel binning/gridding with no extra time-dependent frequency transformations should use mode='cubedata' and the output frame will be 'NOFRAME'. Example: MSs transformed into time-dependent frames, with the desire to image it as is.

  3. If the mstransformed MS is labeled as 'REST', the base frame of the output image will be 'REST', and this is another way to generate images with no extra time-dependent frequency transformations.


3.4.2.4 Image definition by the user (spectral axis)

  1. start : The units of 'start' will decide whether the user means 'channel index', 'frequency' or 'velocity'.

    start=3 : channel 3 of the first spw in the selected list ( irrespective of chan selection in 'spw' )
    start='1.2GHz' : start frequency. channels are equidistant in frequency.
    start='15 km/s' : start velocity.
    $      $ if 'veltype'='RADIO', channels are equidistant in both freq and velo.
    $      $ if 'veltype'='OPTICAL / Z', channels are equidistant in velocity but not in frequency.

  2. step : If 'step' has units, it has to be the same units as start. If specified as an integer, it is taken as N x data_chan_width. For irregular channels in either frequency or velocity, a reasonable default for data_chan_width will be calculated.

  3. frame : Reference frame in which to interpret start : LSRK, LSRD, BARY, GEO, TOPO, GALACTO, LGROUP, CMB This is also the frame to which the conversion layer will be set, at the end of 'clean', for mode='cube'. For mode='cubesrc', the option of specifying start in the SOURCE frame will also be allowed.

  4. veltype : Velocity option in which to interpret start if units are 'km/s' : RADIO, Z, RATIO, BETA, GAMMA, OPTICAL.

  5. restfreq : A vector of rest-frequencies, which will be encoded in the output image. If this parameter is left empty, the (list of) rest-frequencies encoded in the SOURCE table corresponding to the field being imaged, will be used.

    The first rest-frequency in the list will be used to interpret start when its units indicate a velocity specification.

  6. sysvel : Systemic Velocity of a source (only for mode='cubesrc').
  7. sysvelframe : Frequency frame w.r.to which the systemic velocity is specified (only for mode='cubesrc')

  8. nchan : Number of image channels.

  9. interpolation : Relevant when image channel widths > data channel widths, and/or 'start' is offset from the data start. This is irrespective of whether time-dependent frame conversions happen or not. It is a No-Op *only* when 'start' and 'step' are aligned between the data chans and image chans and no time-dependent frequency relabeling is needed.

NOTE : Data selection via 'spw' : The user should select a range larger than what the image will need, and not try to fine-tune the list of channels. The channel mapping and binning process (section 3.4.2) will pick and grid only those data channels that actually map into image channels. We will ensure this is optimized for performance.


3.4.2.5 Using output images from 'clean'

Images from 'clean' will have LSRK or SOURCE or NOFRAME or REST as the base frame. Images with LSRK as the base frame will have a conversion layer set to whatever frame the user specified the input in. The spectral axis of the base frame is always encoded in frequency, in the output images. A regularly spaced set of frequencies is represented by the start/step/nchan listed in the image header. An irregularly spaced set of frequencies is encoded as a tabular axis.

Some details about the conversion layer are

  1. For LSRK output images, a conversion layer will be set to the 'frame' in which the user defined the image channelization. This allows on-the-fly relabeling of the image frequencies from LSRK to 'frame' when read out as world coordinates.

  2. The Viewer will, by default, honor the conversion layer present in an image, and relabel image frequencies on-the-fly while displaying the spectral coordinate. The Viewer also has options to temporarily change the conversion layer to any frequency frame or velocity convention.

  3. Note that conversion layers from LSRK to TOPO/GEO ( time-independent frame to time-dependent frame ) will be tied to one particular time during the observation. Our convention is the start time of the dataset being imaged.

  4. The conversion layer attached to the image can be changed outside the Viewer, using the imreframe task, or a tool script that uses cs.setconversiontype().

  5. Tool level scripts using the imageanalysis and coordinatesystem ia, cs modules can be used to extract lists of frequencies or velocities in any spectral frame and velocity convention. With a conversion layer, csys = ia.coordsys(); csys.toworld( [0,0,0,5] ) will give the frequency of channel 5 in the frame of the conversion layer. With no conversion layer, it will list channel 5 in the base frame of the image ( i.e. LSRK ). Velocities can be read out using csys helper functions : csys.(set)restfrequency(XXX), csys.frequencytovelocity( 1.5, 'GHz', 'RADIO', 'km/s ) . Several other spectral axis relabeling options are possible in combination with the measured me module.

For SOURCE, REST and NOFRAME output images, no conversion layer will be set by clean.

CASA Images can finally be exported to the FITS format, during which, frame conversions are hard-coded. ( Note : This is currently impossible as the exportfits task ignores conversion layers and needs the base frame of the image to be the same as what is desired in the FITS image. CAS-5532 has been made to fix this. ).

Image channels can be regridded using the imregrid task, if the user needs an explicit regrid instead of only frequency-axis relabeling.


3.4.2.6 Frequency types/conversions in casacore::Measures

This section lists the frame definitions and conversion options implemented in the CASACore package ( which CASApy uses ). This information is mostly copied from the code documentation of the MFrequency, MRadialVelocity and MDoppler classes of the Measures module.

Details on how we use this in Imager are also described briefly when relevant.

The different frequency types are: 

    REST -- Rest frequency
    LSRD -- Local Standard of Rest (J2000) -- as the dynamical definition (IAU, [9,12,7] km/s in galactic coordinates)
    LSRK -- LSR as a kinematical (radio) definition -- 20.0 km/s in direction ra,dec = [270,+30] deg (B1900.0)
    BARY -- Barycentric (J2000)
    GEO --- Geocentric
    TOPO -- Topocentric
    GALACTO -- Galacto centric (with rotation of 220 km/s in direction l,b = [90,0] deg.
    LGROUP -- Local group velocity -- 308km/s towards l,b = [105,-7] deg (F. Ghigo)
    CMB -- CMB velocity -- 369.5km/s towards l,b = [264.4, 48.4] deg (F. Ghigo)
    DEFAULT = LSRK

The different types of Doppler (with F = f/f0, the frequency ratio), are: 

    Z = (-1 + 1/F)
    RATIO = (F) *
    RADIO = (1 - F)
    OPTICAL == Z
    BETA = ((1 - F2)/(1 + F2))
    GAMMA = ((1 + F2)/2F) *
    RELATIVISTIC == BETA (== v/c)
    DEFAULT == RADIO
Note that the ones with an '*' have no real interpretation (although the calculation will proceed) if given as a velocity.

An MFrequency can be generated from a simple value (or an MFrequency object), which is then interpreted as a frequency in Hz, and a reference, with an LSRK type as default.

An MRadialVelocity can be generated from a simple value (or an MVRadialVelocity object), which is then interpreted as a RadialVelocity in m/s, and a reference, with an LSRK type as default.

An MDoppler can be generated from a simple value (or an MVDoppler), which is then interpreted as a Doppler ratio, and a reference, with a RADIO type as default. MDoppler can be created from an MFrequency or MRadialVelocity Object, and converted to them too.

All three above objects can also be generated from a Quantity (value/unit pair), where the interpretation depends on the dimensionality of the Quantity: 

----------------------------------------------------------
MFrequency : 
----------------------------------------------------------
    time (e.g. s): period
    frequency (e.g. Hz): frequency
    angular frequency (e.g. arcmin/s): angular frequency
    length (e.g. cm): wavelength
    inverse length (e.g. mm-1): wave number
    energy (e.g. J.s): energy (i.e. h.nu)
    momentum (e.g. kg.m): m.c/h
-----------------------------------------------------------
MVelocity : 
----------------------------------------------------------
    velocity (e.g. AU/a)
----------------------------------------------------------
MDoppler :
----------------------------------------------------------
    None: a Doppler ratio
    Velocity: Doppler ratio calculated by dividing with c
----------------------------------------------------------

Imager allows the user to specify the desired image channelization via start, step, veltype parameters which can be in any of these units. This is internally converted to a Quantity object, from which an MFrequency (frequency in Hz associated with a reference frame) is constructed to populate the spectral axis of the image. The reference frame is the base frame of the image, which depends on the mode of cube imaging (LSRK, SOURCE, REST). In the case of mode='cubedata' and NOFRAME, there is no valid frequency-reference, and channel frequencies are interpreted literally from the start,step,veltype,restfreq parameters.

Conversion between the different types is done with the standard MeasConvert class (MFrequency::Convert, MRadialVelocity::Convert, MDoppler::Convert). This is what is encoded in the conversion layer of CASA Images.

Some conversions are only possible if the following frame information is available.

  1. Conversion to/from REST needs Radial Velocity information. The sysvel parameter in mode='cubesrc' will be used for this. For an MS already at REST, no conversions are needed.
  2. Conversion to/from TOPO and GEO needs Epoch information. This is set in the conversion layer for mode='cube' as the start time of the MS and can be modified via the imreframe task with outframe='TOPO' or 'GEO' and subparameter epoch.
  3. Conversion to/from TOPO needs Position information. This is read from the input MS, or Image header.
  4. All conversions need Direction information. This is the image center from the Image header.

Some Misc. info :

Caution: Conversion between the different frequencies can, due to relativistic effects, only be done approximately for very high (order c) radial velocities (shifted frequencies); A better approach would be to start from radial velocities and a rest frequency;

Caution: For large radial velocities (of order c) the conversions are not precise, and not completely reversable, due to unknown transverse velocities, and the additive way in which corrections are done; They are correct to first order wrt relativistic effects

If the Doppler shift is known (e.g. from another spectral line), the REST frequency can be determined with the toREST() member. Dopplers do not need a reference frame.

next up previous contents
Next: 4. CASA Imager Memory Up: Imaging Algorithms in CASA Previous: 2. Currently available Imaging   Contents
R. V. Urvashi 2013-10-22