This page gives an overview of how to do operations on images using the CASA tasks immath, imstat, and imdev.

Mathematical Operations on an Image (immath)

The inputs are:

In all cases, outfile must be supplied with the name of the new output file to create. The mode parameter selects what immath is to do. The default mode='evalexpr' lets the user specify a mathematical operation to carry out on one or more input images. The sub-parameter expr contains the Lattice Expression Language (LEL) string describing the image operations based on the images in the imagename parameter.

Mask specification is done using the mask parameter. This can optionally contain an on-the-fly mask expression (in LEL) or point to an image with a pixel mask. See above for more on the use of the mask parameter. See also the page on LEL strings. Sometimes, one would like to use a flat image (e.g. a moment image) mask to be applied to an entire cube. The stretch=True subparameter in mask allows one to expand the mask to all planes of the cube.

Region selection is carried out through the region and box parameters (see above). Image plane selection is controlled by chans and stokes parameters. For mode='evalexpr', the standard usage for specifying images to be used in the LEL expression is to provide them as a list in the imagename parameter, and then access there in the LEL expression by the names IM0, IM1, .... E.g.,

would subtract the second image given from the first. For the special modes 'spix', 'pola', 'poli', the required images for the given operation are to be provided in imagename (sometimes in a particular order).

The parameter imagemd can be set to a image name from which the header meta-information is to be used for the output image. If left unset the task may pick any of the input image headers, so it is better to define this parameter. In fact, the image specified in imagemd can be any image, even an image that is not part of the calculations in immath.

Detailed examples are given below.

Examples for immath

In the following, we show a few examples of immath. Note that the image names in the expr are assumed to refer to existing image files in the current working directory.

 

Simple math

Select a single plane (channel 22) of the 3-D cube:

Double all values in our image:

Square all values in our image:

You can do other mathematical operations on an image (e.g. trigonometric functions), as well as use scalars results from an image (e.g. max, min, median, mean, variance). You also have access to constants such as e() and pi() (which are doubles internally, while most images are floats). For example: Take the sine of an image:

Note again that the units are again kept as they were. Select a single plane (channel 22) of the 3-D cube and subtract it from the original image:

Note that in this example the 2-D plane gets extended in the third dimension and the 2-D values are applied to each plane in the 3-D cube. Select and save the inner 1/4 of an image for channels 40,42,44 as well as channels 10 and below:

If chans selects more than one channel then the output image has a number of channels given by the span from the lowest and highest channel selected in chans. In the example above, it will have 45 channels. The ones not selected will be masked in the output cube. If we had set chans='40,42,44' then there would be 5 output channels corresponding to channels 40,41,42,43,44 of the MS with 41,43 masked. Also, the chans='<10' selects channels 0–9. Note that the chans syntax allows the operators '<', '<=', '>', '>'. For example,

do the same thing. Divide an image by another, with a threshold on one of the images:

Polarization manipulation

It is helpful to extract the Stokes planes from the cube into individual images:

Extract linearly polarized intensity and polarization position angle images:

Create a fractional linear polarization image:

Create a polarized intensity image:

Toolkit Tricks: The following uses the toolkit. You can make a complex linear polarization (Q+iU) image using the imagepol tool:

You can now display this in the viewer, in particular overlay this over the intensity raster with the intensity contours. When you load the image, use the LEL:

which is entered into the LEL box at the bottom of the Load Data menu.

 

Using masks in immath

The mask parameter is used inside immath to apply a mask to all the images used in expr before calculations are done (if you are curious, it uses the ia.subimage tool method to make virtual images that are then input in the LEL to the ia.imagecalc method).

For example, let’s assume that we have made a single channel image using clean:

There is now a file ngc5921.demo.chan22.cleanimg.mask that is an image with values 1.0 inside the cleanbox region and 0.0 outside. We can use this to mask the clean image:

Toolkit Tricks: Note that there are also pixel masks that can be contained in each image. These are Boolean masks, and are implicitly used in the calculation for each image in expr. If you want to use the mask in a different image not in expr, try it in mask:

Note that nominally the axes of the mask must be congruent to the axes of the images in expr. However, one exception is that the image in mask can have fewer axes (but not axes that exist but are of the wrong lengths). In this case, immath will extend the missing axes to cover the range in the images in expr. Thus, you can apply a mask made from a single channel to a whole cube.

 

Computing image statistics (imstat)

The imstat task will calculate statistics on a region of an image, and return the results as a return value in a Python dictionary. The inputs are:

Area selection can be done using region and mask parameters. Plane selection is controlled by chans and stokes. The parameter axes will select the dimensions that the statistics is calculated over. Typical data cubes have axes like: RA axis 0, DEC axis 1, Velocity axis 2. So, e.g. axes=[0,1] would be the most common setting to calculate statistics per spectral channel.

A typical output of imstat on a cube with axes=[0,1] and algorithm='classic' (default) looks like:

where the header information provides the specifications of the data that were selected followed by the table with the frequency values of the lanes, the plane numbers, Npts (the number of pixels per plane), and the Sum, Median, RMS, Standard deviations, Minimum, and Maximum of the pixel values for each plane. Similar output is provided when the data is averaged over different axes. The logger output can also be written into or appended to a log file for further processing elsewhere (logfile parameter).

imstat has access to different statistics algorithms. Most of them represent different ways on how to treat distributions that are not Gaussian, in particular to eliminate outlier values from the statistics. Available algorithms are CLASSIC, where all unmasked pixels are used, FIT-HALF, where one (good) half of the distribution is being mirrored across a central value, HINGES-FENCES, where the inner quartiles plus a ’fence’ data portion is being used, and CHAUVENET, which includes values based on the number of standard deviations from the mean. For more information, see the inline help of the imstat task.

Using the task return value

The contents of the return value of imstat are in a Python dictionary of key-value sets. For example,

will assign this to the Python variable xstat. The keys for xstat are outlined on the imstat page.

For example, an imstat call might be

In the terminal window, imstat reports:

The return value in xstat is

For example, to extract the standard deviation as a number

Examples for imstat

To extract statistics for an image:

In a box around the brightest component:

 

Computing a deviation image (imdev)

The imdev task produces an output image whose value in each pixel represents the "error" or "deviation" in the input image at the corresponding pixel. The output image has the same dimensions and coordinate system as the input image, or the selected region of the input image.

Area selection can be done using region and mask parameters. Plane selection is controlled by chans and stokes. Statistics are computed spatially: a deviation image is computed independently for each channel/Stokes plane. If the outfile parameter is left blank, the task returns an image tool referencing the resulting image; otherwise the resulting image is written to disk.

The statistic to be computed is selected using the stattype parameter. Allowed statistics are:

The chosen statistic is calculated around a set of grid points (pixels) across the input image. The grid parameter is an array of two positive integers [x,y] specifying the grid spacing in the x and y direction. The anchor parameter specifies a pixel used to "anchor" the grid; the default is to use the reference pixel of the input image. Using the xlength and ylength parameters, the user may specify the dimensions of a rectangular or circular region around each grid point that defines which surrounding pixels are used in the statistic computation for individual grid points. If the ylength parameter is an empty string, then a circle centered on each grid point with diameter provided by xlength is used for the statistic computation. Otherwise, a rectangular region centered on each grid point with dimensions of xlength x ylength is used. These two parameters may be specified as valid quantities with recognized units (e.g., "4arcsec"). They may also be specified as numerical values, in which case the unit is assumed to be pixels (e.g., an input value of 4 is equivalent to an input value of "4pix").

The chosen statistic is calculated at every grid point in the input image, and the result is reflected at the corresponding pixel of the output image. Values at all other pixels in the output image are determined by interpolating across the grid points. The interp parameter can be used to specify the algorithm for the computation. Available algorithms are CLASSIC, where all unmasked pixels are used, and CHAUVENET, which includes values based on the number of standard deviations from the mean.

Examples for imdev

Compute a "standard deviation" image using grid-spacing of 4 arcsec in the X direction and 5 arcsec in the Y direction, with linear interpolation to compute values at non-grid-pixels. Compute the standard deviation in a box of 20 x 25 arcsec.

Compute an image showing median absolute deviation (MAD) across the image, with MAD converted to an equivalent RMS value. Anchor the grid at a specific pixel [1000,1000] with grid-spacing of 10 pixels, and use circles of diameter 30 pixels for the statistical computation. Calculate the statistic using the z-score/Chauvenet algorithm by fixing the maximum z-score to determine outliers to 5. Use cubic interpolation to determine the value at non-grid-pixels. Have the task return a pointer to the output image.