Math Operations / Statistics

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:

#  immath :: Perform math operations on images
imagename           =         ''        #  a list of input images
mode                = 'evalexpr'        #  mode for math operation (evalexpr, spix, pola, poli)
     expr           =         ''        #  Mathematical expression using images
     varnames       =         ''        #  a list of variable names to use with the image files

outfile             = '' #  File where the output is saved
mask                =         ''        #  Mask to use. Default is none.
region              =         ''        #  Region selection. 
                                        #   Default is to use the full image.
box                 =         ''        #  Rectangular region to
                                        #   select in direction plane.
                                        #    Default is to use the
                                        #   entire direction plane.
chans               =         ''        #  Channels to use. 
                                        #   Default is to use all channels.
stokes              =         ''        #  Stokes planes to use. 
                                        #   Default is to use all Stokes planes.
imagemd             =         ''        #  An image name from which metadata should be copied. The input
                                        #   can be either an image listed under imagename or any other
                                        #   image on disk. Leaving this parameter unset may copy header
                                        #   metadata from any of the input images, which
                                        #   one is not guaranteed.

Alert: immath does not convert any brightness units, e.g. from Jy/beam to K or vice versa. The user is responsible for making sure the images are consistent with the values in the header and the image. It is not advisable to mix input images that are in different units or have different beam sizes.

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

ALERT: For mode=’pola’ you MUST call as a function as in the example below, giving the parameters as arguments, or immath will fail.

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:

       outfile='ngc5921.demo.chan22double.image' )

Square all values in our image:

       outfile='ngc5921.demo.chan22squared.image' )

NOTE: The units in the output image are still claimed to be “Jy/beam”, i.e. immath will not correctly scale the units in the image for non-linear cases like this. Beware!

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,

chans = '<17,>79'
chans = '<=16,>=80'

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:

imagename = '3C129BC.clean.image'
outfile='3C129BC.I'; expr='IM0'; stokes='I'; immath();
outfile='3C129BC.Q'; expr='IM0'; stokes='Q'; immath();
outfile='3C129BC.U'; expr='IM0'; stokes='U'; immath();
outfile='3C129BC.V'; expr='IM0'; stokes='V'; immath();

Extract linearly polarized intensity and polarization position angle images:

immath(stokes='', outfile='3C129BC.P', mode='poli',
       imagename=['3C129BC.Q','3C129BC.U'], sigma='0.0mJy/beam');
immath(stokes='', outfile='3C129BC.X', mode='pola',
       imagename=['3C129BC.Q','3C129BC.U'], sigma='0.0mJy/beam');

ALERT: For mode=’pola’ you MUST call as a function as in this example (giving the parameters as arguments) or immath will fail.

Create a fractional linear polarization image:

default( 'immath')
imagename = ['3C129BC.I','3C129BC.Q','3C129BC.U']
expr='sqrt((IM1^2 + IM2^2)/IM0^2)'

Create a polarized intensity image:

default( 'immath')
imagename = ['3C129BC.Q','3C129BC.U','3C129BC.V']
expr='sqrt(IM0^2 + IM1^2 + IM2^2)'

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

  # Make an imagepol tool and open the clean image
  potool = casac.homefinder.find_home_by_name('imagepolHome')
  po = potool.create()'3C129BC.clean.image')
  # Use complexlinpol to make a Q+iU image

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:

vis = ''
imagename = 'ngc5921.demo.chan22.cleanimg'
mode = 'channel'
nchan = 1
start = 22
step = 1
field = ''
spw = ''
imsize = [256,256]
cell = [15.,15.]
psfalg = 'clark'
gain = 0.1
niter = 6000
weighting = 'briggs'
rmode = 'norm'
robust = 0.5

mask = [108,108,148,148]

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:

imagename = 'ngc5921.demo.chan22.cleanimg.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:

# First make a pixel mask inside ngc5921.demo.chan22.cleanimg.mask'ngc5921.demo.chan22.cleanimg.mask')
# There is now a 'mask0' mask in this image as reported by the summary

# Now apply this pixel mask in immath

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.

# drop degenerate stokes and freq axes from mask image'ngc5921.demo.chan22.cleanimg.mask')
im2 = ia.subimage(outfile='ngc5921.demo.chan22.cleanimg.mymask',dropdeg=True)
# mymask has only RA and Dec axes

# Now apply this mask to the 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:

#  imstat :: Displays statistical information from an image or image region
imagename           =         ''        #  Name of the input image
axes                =         -1        #  List of axes to evaluate statistics over. Default is
                                        #   all axes.
region              =         ''        #  Image Region or name. Use Viewer
box                 =         ''        #  Select one or more box regions
chans               =         ''        #  Select the channel(spectral) range. 
stokes              =         ''        #  Stokes params to image (I,IV,IQU,IQUV). Default '' =>
                                        #   include all
listit              =       True        #  Print stats and bounding box to logger?
verbose             =      False        #  Print additional messages to logger?
mask                =         ''        #  Mask to use. Default is none.
logfile             =         ''        #  Name of file to write fit results.
algorithm           =  'classic'        #  Algorithm to use. Supported values are 'chauvenet',
                                        #   'classic', 'fit-half', and 'hinges-fences'. Minimum
                                        #   match is supported.
     clmethod       =     'auto'        #  Method to use for calculating classical statistics.
                                        #   Supported methods are 'auto', 'tiled', and
                                        #   'framework'. Ignored if algorithm is not 'classic'.

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:

No region specified. Using full positional plane.
Using all spectral channels.
Using polarizations ALL
Determining stats for image IRC10216_HC3N.cube_r0.5.image
Set region from supplied region record
Statistics calculated using Classic algorithm
Regions ---
         -- bottom-left corner (pixel) [blc]:  [0, 0, 0, 0]
         -- top-right corner (pixel) [trc]:    [299, 299, 0, 63]
         -- bottom-left corner (world) [blcf]: 09:48:01.492, +, I, 3.63994e+10Hz
         -- top-right corner (world) [trcf]:   09:47:53.299, +, I, 3.63915e+10Hz
No region specified. Using full positional plane.
Using all spectral channels.
Using polarizations ALL
Selected bounding box :
    [0, 0, 0, 0] to [299, 299, 0, 63]  (09:48:01.492, +, I, 3.63994e+10Hz to 09:47:53.299, +, I, 3.63915e+10Hz)
#        Frequency  Frequency(Plane) Npts          Sum           Mean          Rms           Std dev       Minimum       Maximum     
  3.63993552e+10                  0  9.000000e+04  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
  3.63992302e+10                  1  9.000000e+04  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
  3.63991052e+10                  2  9.000000e+04  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
  3.63989802e+10                  3  9.000000e+04  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
  3.63988551e+10                  4  9.000000e+04  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
  3.63987301e+10                  5  9.000000e+04  6.069948e-01  6.744386e-06  1.534640e-03  1.534634e-03 -6.355108e-03  6.166496e-03
  3.63986051e+10                  6  9.000000e+04  2.711720e-01  3.013023e-06  1.538071e-03  1.538077e-03 -6.165663e-03  5.862981e-03
  3.63984801e+10                  7  9.000000e+04  2.501259e-01  2.779177e-06  1.578049e-03  1.578056e-03 -6.771976e-03  6.272645e-03
  3.63983551e+10                  8  9.000000e+04 -3.706732e-01 -4.118591e-06  1.607191e-03  1.607194e-03 -8.871284e-03  6.591001e-03

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,

xstat = imstat()

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

 imagename = 'ngc5921.demo.cleanimg.image'  #  The NGC5921 image cube
 box       = '108,108,148,148'              #  20 pixels around the center
 chans     = '21'                           #  channel 21

 xstat = imstat()

In the terminal window, imstat reports:

Statistics on  ngc5921.usecase.clean.image

Region ---
   -- bottom-left corner (pixel) [blc]: [108, 108, 0, 21]
   -- top-right corner (pixel) [trc]:   [148, 148, 0, 21]
   -- bottom-left corner (world) [blcf]: 15:22:20.076, +, I, 1.41332e+09Hz
   -- top-right corner( world) [trcf]: 15:21:39.919, +, I, 1.41332e+09Hz

Values --
   -- flux [flux]:              0.111799236126
   -- number of points [npts]:  1681.0
   -- maximum value [max]:      0.029451508075
   -- minimum value [min]:     -0.00612453464419
   -- position of max value (pixel) [maxpos]:  [124, 131, 0, 21]
   -- position of min value (pixel) [minpos]:  [142, 110, 0, 21]
   -- position of max value (world) [maxposf]: 15:22:04.016, +, I, 1.41332e+09Hz
   -- position of min value (world) [minposf]: 15:21:45.947, +, I, 1.41332e+09Hz
   -- Sum of pixel values [sum]: 1.32267159822
   -- Sum of squared pixel values [sumsq]: 0.0284534543692
Statistics ---
   -- Mean of the pixel values [mean]:       0.000786836167885
   -- Standard deviation of the Mean [sigma]: 0.00403944306904
   -- Root mean square [rms]:               0.00411418313161
   -- Median of the pixel values [median]:     0.000137259965413
   -- Median of the deviations [medabsdevmed]:       0.00152346317191
   -- Quartile [quartile]:                       0.00305395200849

The return value in xstat is

CASA <152>: xstat
{'blc': array([108, 108,   0,  21]),
 'blcf': '15:22:20.076, +, I, 1.41332e+09Hz',
 'flux': array([ 0.11179924]),
 'max': array([ 0.02945151]),
 'maxpos': array([124, 131,   0,  21]),
 'maxposf': '15:22:04.016, +, I, 1.41332e+09Hz',
 'mean': array([ 0.00078684]),
 'medabsdevmed': array([ 0.00152346]),
 'median': array([ 0.00013726]),
 'min': array([-0.00612453]),
 'minpos': array([142, 110,   0,  21]),
 'minposf': '15:21:45.947, +, I, 1.41332e+09Hz',
 'npts': array([ 1681.]),
 'quartile': array([ 0.00305395]),
 'rms': array([ 0.00411418]),
 'sigma': array([ 0.00403944]),
 'sum': array([ 1.3226716]),
 'sumsq': array([ 0.02845345]),
 'trc': array([148, 148,   0,  21]),
 'trcf': '15:21:39.919, +, I, 1.41332e+09Hz'}

ALERT: The return dictionary currently includes NumPy array values, which have to be accessed by an array index to get the array value. To access these dictionary elements, use the standard Python dictionary syntax, e.g. xstat[][]

For example, to extract the standard deviation as a number

mystddev = xstat['sigma'][0]
print 'Sigma = '+str(xstat['sigma'][0])

Examples for imstat

To extract statistics for an image:

xstat = imstat('b1608.demo.clean2.image')
# Printing out some of these
  print 'Max   = '+str(xstat['max'][0])
  print 'Sigma = '+str(xstat['sigma'][0])
# results:
# Max   = 0.016796965152
# Sigma = 0.00033631979385

In a box around the brightest component:

xstat_A = imstat('b1608.demo.clean2.image',box='124,125,132,133')
# Printing out some of these
  print 'Comp A Max Flux = '+str(xstat_A['max'][0])
  print 'Comp A Max X,Y  = ('+str(xstat_A['maxpos'][0])+','+str(xstat_A['maxpos'][1])+')'
# results:
# Comp A Max Flux = 0.016796965152
# Comp A Max X,Y  = (128,129)


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 as the selected region of the input image.

# imdev :: Create an image that can represent the statistical deviations of the input image.
imagename          =          ''        # Input image name
outfile            =          ''        # Output image file name. If left blank (the default), no image is written but a new image
                                        # tool referencing the collapsed image is returned.
region             =          ''        # Region selection. Default is to use the full image.
box                =          ''        # Rectangular region(s) to select in direction plane. Default is to use the entire
                                        # direction plane.
chans              =          ''        # Channels to use. Default is to use all channels.
stokes             =          ''        # Stokes planes to use. Default is to use all Stokes planes.
mask               =          ''        # Mask to use. Default setting is none.
overwrite          =       False        # Overwrite (unprompted) pre-existing output file? Ignored if "outfile" is left blank.
grid               =      [1, 1]        # x,y grid spacing. Array of exactly two positive integers.
anchor             =       'ref'        # x,y anchor pixel location. Either "ref" to use the image reference pixel, or an array of
                                        # exactly two integers.
xlength            =      '1pix'        # Either x coordinate length of box, or diameter of circle. Circle is used if ylength is
                                        # empty string.
ylength            =      '1pix'        # y coordinate length of box. Use a circle if ylength is empty string.
interp             =     'cubic'        # Interpolation algorithm to use. One of "nearest", "linear", "cubic", or "lanczos".
                                        # Minimum match supported.
stattype           =     'sigma'        # Statistic to compute. See full description for supported statistics.
statalg            =   'classic'        # Statistics computation algorithm to use. Supported values are "chauvenet" and "classic",
                                        # Minimum match is supported.

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:

iqr                      inner quartile range (q3 - q1)
max                      maximum
mean                     mean
medabsdevmed, madm       median absolute deviation from the median
median                   median
min                      minimum
npts                     number of points
q1                       first quartile
q3                       third quartile
rms                      rms
sigma, std               standard deviation
sumsq                    sum of squares
sum                      sum
var                      variance
xmadm                    median absolute deviation from the median multipied by x, where x is the reciprocal of Phi^-1(3/4),
                         where Phi^-1 is the reciprocal of the quantile function. Numerically, x = 1.482602218505602. See, eg,

The chosen statistic is calculated around a set of grid points (pixels) across the input image with grid spacing specified by the grid parameter. The size and shape of the region used to compute the statistic at each grid point is specified by the xlength and ylength parameters. If ylength is an empty string, then the region used is a circle centered on each grid point with diameter provided by xlength. Otherwise, a rectangular region with dimensions of xlength by ylength is used. These two parameters may be specified as valid quantities with recognized units (e.g., "4arcsec" or "4pix"). They may also be specified as numerical values, in which case the unit is assumed to be pixels.

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 using the interpolation scheme given by the input parameter interp. The statalg parameter specifies the algorithm for the statistics 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.

imdev("my.image", "std.image", grid=[4,5], xlength="20arcsec", ylength="25arcsec", stattype="sigma", interp="linear", statalg="classic")

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.

myim = imdev("my.image", anchor=[1000,1000], grid=[10,10], xlength=30, ylength='', stattype="xmadm", interp="cubic", statalg="chauvenet", zscore=5)