Perform arithmetic on images
dmimgcalc infile infile2 outfile operation [weight] [weight2] [lookupTab] [clobber] [verbose]
`dmimgcalc' is used to combine multiple images using a wide range of mathematical operations. Rather than being restricted to simple operations, (almost) arbitrarily complex expressions can be used. The tool can also be used to compare two images on a pixel-by-pixel basis. The input images must have the same sizes and are combined using logical coordinates.
Simple combinations of two images
Using the operation, weight, and weight2 parameters, one or two images can be combined - using addition, subtraction, multiplication, or division - with optional weights. The output image is defined as
image1*weight <operation> image2*weight2
when <operation> is one of "add", "sub", "mul", or "div". The data type of the output image will be set to double if weighting is used (either weight or weight2 is not equal to 1), one integer image is dividing another integer image, or the two images have different data types, otherwise the data type will match that of the input images.
Comparing two images
If operation is set to "tst" then the two images are checked for equality - i.e. that matching pixels in the two images are the same - and the exit status value set accordingly (0 if equal, 1 if not). Note that in this case there is no output image so outfile can be set to "none".
Adding or subtracting a constant
To add or subtract a constant from an image, there must be a single input file (not a stack) and infile2=none. Then set the operation as shown:
# addition unix% dmimgcalc infile=img.fits op="imgout=img1+10" \ out="img_add.fits" mode=h # subtraction unix% dmimgcalc infile=img.fits op="imgout=img1-2" \ out="img_sub.fits" mode=h
Multiplying an image by a constant
There are two different methods for multiplying an image by a constant. In both cases, there must be a single input file (not a stack) and infile2=none.
- Set "weight" to the multiplication factor and the operation to any of "add", "sub", "mul", or "div". The output image is the input image multiplied by the weight.
- Set the operation parameter as "img1*N" where N is the multiplication factor.
These two multiplication commands yield comparable results:
unix% dmimgcalc infile=img.fits infile2=none op=add \ out=img_mult_2.fits weight=5 unix% dmimgcalc infile=img.fits op="imgout=img1*5" \ out="img_mult.fits" mode=h
General arithmetic combinations of many images
The previous syntax is adequate for simple image combinations. More complex operations - such as combining more than two images, converting the data type of an image, or applying mathematical functions to images - can also be specified. In this case the list of input images is given as a stack to the infile parameter (see "ahelp stack") and the individual images are referenced by the labels "img<n>", where "img1" refers to the first image and "imgn" the nth image in the stack.
In this mode the operation parameter is set to
where the allowed syntax for everything to the right of "imgout=" is described in "ahelp dmmath". As an example, setting
will set the output image to be the sum of the square roots of the first two images divided by the third image. Further examples are given in the Examples section below. Note that as of CIAO 4.15, the operation parameter can be set to a stack - such as op=@coms - where the file in question (e.g. coms) must start with the string "imgout=" (without the quotes).
In this mode the weight and weight2 parameters are ignored.
The tool will operate on corresponding pixels from the two images regardless of any (conflicting) WCS information in the header. Thus the user must ensure that the two images must be aligned in logical array coordinates prior to using the tool. If infile has WCS information, that information will be propagated to the output file. If the WCS in infile does not match the WCS in infile2, a warning will be printed.
dmimgcalc acis.fits exposure.fits normalized.fits div
The file normalized.fits is set to the results of dividing acis.fits by exposure.fits.
dmimgcalc acis1.fits acis2.fits output.fits add
Add images acis1.fits and acis2.fits and store in output.fits.
dmimgcalc acis1.fits acis2.fits output.fits add weight=2.5 weight2=3
Add acis1.fits*2.5 to acis2.fits*3 and store in output.fits.
dmimgcalc acis1.fits none output.fits add weight=3.2
Multiply acis1 by 3.2 and store in output.fits. Note that the value for the operation parameter is ignored (as long as it is set to a legal value).
dmimgcalc can also be used for more complicated image combinations, as shown in the following examples.
dmimgcalc a.fits none sa.fits op="imgout=sqrt(fabs(img1))"
Here we use the ability to specify aribtrary expressions to create sa.fits which contains the square root of the absolute pixel values in the input image (a.fits).
Since the value of the operation parameter begins with 'imgout=' we have to include the parameter name - otherwise the tool would think you were setting the parameter imgout to a value which would result in an error.
dmimgcalc img.fits none dimg.fits op="imgout=(double)img1"
Here we use dmimgcalc to convert the input image (img.fits) into double format (dimg.fits).
dmimgcalc real8.fits none real4.fits op="imgout=(float)img1"
Here dmimgcalc is used to convert the input image from double format to float.
dmimgcalc a.fits,b.fits,c.fits none sum.fits op="imgout=img1+img2+img3"
Here we use dmimgcalc to add the three images together to form sum.fits.
dmimgcalc A,B,C,D,..,Z none mean.fits op="imgout=((img1+img2+..+img26)/26.0)"
Average a stack of 26 images. For another method of averaging images, see "ahelp dmimgfilt".
dmimgcalc "a.fits,b.fits" none combined.fits op="imgout=(((img1*img1_exposure)+(img2*img2_exposure))/(img1_exposure+i mg2_exposure))"
Here we combine a stack of input files (a.fits and b.fits) based on the expression given to the operation parameter to create combined.fits. The expression used is the weighted sum of two images, using the exposure times of the two images as the weights.
In this mode, the stack of files input as the infile parameter (so a.fits and b.fits in the example) are referred to be "img<n>" in the expression (<n> starts from 1, so "img1" refers to a.fits in the example). The output image is referenced by "imgout", and the "img1_exposure" terms are replaced by the value of the "EXPOSURE" keyword from the header of the images (in this case a.fits). So, if the EXPOSURE values of a.fits and b.fits were 12012.45 and 24034.9 respectively, the output image (combined.fits) would be equal to
( (a.fits*12012.45)+(b.fits*24034.9) ) / (12012.45+24034.9)
cat << EOM > expression.lis imgout=(\ ((img1*img1_exposure)+(img2*img2_exposure))/\ (img1_exposure+img2_exposure)) EOM dmimgcalc "a.fits,b.fits" none combined.fits email@example.com
This is the same as the previous example, except now the long expression has been saved into an external ASCII file and is being read in using the CIAO stack syntax. It is important to note how each of the first two lines ends with a trailing back-slash, "\". This "continuation character" means that the next line is meant to continue the current line. Using this syntax users can create much longer expressions than can currently be input via the parameter directly. Expressions up to 32,000 characters can be used this way.
dmimgcalc acis1.fits acis2.fits none tst
If the two images (here acis1.fits and acis2.fits) are equal - in that corresponding pixels have the same value in the two images - then the status/exit value is set to 0, otherwise it is set to 1. The way to check this value depends on what shell you are using; two common ways are:
|echo $status||for csh/tcsh users|
|echo $?||for bash/sh users|
Detailed Parameter Descriptions
Parameter=infile (file required filetype=input stacks=yes)
The input file (or files)
If a stack of files is used then the operation parameter should be of the form "imgout=...".
Parameter=infile2 (file required filetype=input)
The second input virtual file specification, or "none". If the former, image must be same size as infile.
Parameter=outfile (file required filetype=output)
The output file name
Parameter=operation (string required stacks=yes)
Arithmetic operation to be performed (add/sub/mul/div/tst/expression).
This parameter can be one of the following:
Valid arithmetic operations
|sub||subtract infile2 image from infile image|
|div||divide infile image by infile2 image|
|tst||if each pixel in infile2 equals the same pixel in infile then set the status flag to 0, otherwise set it to 1. The outfile parameter can be set to "none".|
|imgout=.....||Used to specify an arbitrary expression combining a number of files. The infile2 parameter should be set to "none" when using this form, either directly or via a stack file, and the values of the weight and weight2 parameters are ignored.|
|@file (stack)||As of CIAO 4.15, the arbitrary expression can be read in from a CIAO stack file, such as @expression.lis, where the ASCII file expression.lis starts with "imgout=" (without quotes). This file can contain much-longer expressions than can be used on the command line. To make them readable, a multi-line file can be used as long as the last character of each but the last line ends with a back-slash character (\) to indicate the lines should be concatenated together.|
Parameter=weight (real default=1)
Value by which to multiply infile before operation.
This value is not used when the operation parameter starts with "imgout=".
Parameter=weight2 (real default=1)
Value by which to multiply infile2 before operation
This value is not used when the operation parameter starts with "imgout=".
Parameter=lookupTab (file filetype=input)
Lookup table for header merging
The look-up table defines the rules used to combine header keywords from multiple files, as described in "ahelp merging_rules". Set to "none" to stop the header files from being merged (so the output file will have a copy of the first input file).
Parameter=clobber (boolean default=no)
Overwrite existing file with same output filename?
Parameter=verbose (integer default=0 min=0 max=5)
Determines the level of descriptive output as the tool runs.
Set to "0" for none, up to "5" for maximum description.
Handling Invalid Data
If a division by zero occurs then the output value is set to 0 or NaN, depending on the operation. Attempts to perform any operation on a pixel value of NaN are also set to 0. If these changes occur and the operation was a "simple" one - namely "add", "sub", "mul", or "div" - then a summary message will be displayed at the end of the program.
When images are combined the header keywords from the different files to be merged together following the rules given in the file pointed to in the lookupTab parameter. The rules for the merging - and the syntax for this file - are given in "ahelp merging_rules".
One common occurrence - when using the default set of rules - is for the DETNAM keyword to be changed to the value "Merged", which may cause issues when processing the resulting image with other tools. This tends to happen when combining an observation with an exposure map since the data will likely have a DETNAM field that is set to the value from the event file - e.g. "ACIS-012367" - whilst the exposure map just has the individual chipe name - e.g. "ACIS-7". Since the two values are different the output is set to "Merged". You can set the lookupTab parameter to "NONE" to avoid this merging - in which case the header is taken from the first input file, use a modified copy of the table used to define the merging rules, or use dmhedit to manually set the header keyword to the desired value.
Changes in CIAO 4.15
The operation parameter can now be set to a stack, which is useful when combining many files (where the length of the argument could exceed the 1024 character length limit of the parameter interface). The tool will now error out if given an image file it does not support.
The op parameter can now take a stack list for input. Since stacks support line continuation via a trailing "\" character, they can support arbitrary long expressions. For example
unix% cat my_op.lis imgout=\ (img1*img1_expousre+\ img2*img2_exposure+\ ... img30*img30_exposure)/\ (img1_exposure+\ img2_exposure+\ ... img30_exposure) unix% dmimgcalc @myimgs.lis op=@my_op.lis out=my_out.fits mode=h clob+
Improve error checking for unsupported output data types.
- The operation string cannot exceed 1023 characters.
The operation string is limited to 1023 characters (the imgout= string counts against that 1023 limit). Values longer than 1023 characters cause a string buffer overflow (ie memory corruption) leading to unpredicatable behavior.
- log() of integer values
Using the log or log10 operator on integer pixel values returns an integer value. The value is truncated.
If unsure of the datatype, the value can be explicitly cast to a real-valued number
- Boolean operators - such as < and > - are not working on arrays (e.g. images).
Note that threshold conditions may be more easily applied by using dmimgthresh instead of this tool.
- dmimgcalc corrupts the file history (07 Mar 2008)
dmimgcalc corrupts the history records in the file header, duplicating previous history entries.
- The tool segmentation faults if the input file is a table (28 Feb 2012)
unix% dmimgcalc infile=img.a,img.b infile2= \ outfile=sub.fits op='imgout=img1-img2' warning: ASPTYPE has different value...Merged... BTIMDRFT not present in all input files...FAIL... BTIMNULL not present in all input files...FAIL... BTIMRATE not present in all input files...FAIL... warning: DATAMODE has different value...Merged... warning: DS_IDENT has different value...Merged... omit - FOC_LEN values different more than 1.000000 warning: OBJECT has different value...Merged... warning: OBSERVER has different value...Merged... omit - ROLL_NOM values different more than 1.000000 warning: SEQ_NUM has different value...Merged... warning: TITLE has different value...Merged... # 49970: Received error signal SIGSEGV-segmentation violation. # 49970: An invalid memory reference was made. # 49970: segmentation fault: dmimgcalc (1) is: exit_upon_error->NULL unix% dmlist img.a blocks -------------------------------------------------------------------------------- Dataset: img.a -------------------------------------------------------------------------------- Block Name Type Dimensions -------------------------------------------------------------------------------- Block 1: PRIMARY Null Block 2: EVENTS Table 15 cols x 634221 rows Block 3: GTI7 Table 2 cols x 1 rows
- Mixing scalar and array values (07 Mar 2008)
Users need to be careful when trying to mix scalar and array operations to get the correct results. For example
unix% dmimgcalc img.fits none out.fits op="imgout=5"
Does not produce an image the same size as img.fits whose values are all 5. Rather one gets a 1x1 image with the value=5. To get an image with all values set to 5 we can do something like
unix% dmimgcalc img.fits none out.fits op="imgout=(img1*0)+5" unix% dmimgcalc img.fits none out.fits op="imgout=(img1/img1)*5"
- addresp, dmappend, dmcontour, dmellipse, dmfilth, dmgti, dmimg2jpg, dmimgadapt, dmimgblob, dmimgcalc, dmimgdist, dmimgfilt, dmimghist, dmimghull, dmimglasso, dmimgpick, dmimgpm, dmimgproject, dmimgreproject, dmimgthresh, dmmaskbin, dmmaskfill, dmmath, dmmerge, dmnautilus, dmregrid, dmregrid2, dmstat, dmtcalc, evalpos, hrc_dtfstats, imgmoment, mean_energy_map, mtl_build_gti, pileup_map, reproject_image, reproject_image_grid