GNU Astronomy Utilities

6.2.2 Arithmetic operators

The recognized operators in Arithmetic are listed below. See Reverse polish notation for more on how the operators and operands should be ordered on the command-line. The operands to all operators can be a data array (for example a FITS image) or a number, the output will be an array or number according to the inputs. For example a number multiplied by an array will produce an array. The conditional operators will return pixel, or numerical values of 0 (false) or 1 (true) and stored in an unsigned char data type (see Numeric data types).

+

Addition, so “4 5 +” is equivalent to \(4+5\).

-

Subtraction, so “4 5 -” is equivalent to \(4-5\).

x

Multiplication, so “4 5 x” is equivalent to \(4\times5\).

/

Division, so “4 5 /” is equivalent to \(4/5\).

%

Modulo (remainder), so “3 2 %” is equivalent to \(1\). Note that the modulo operator only works on integer types.

abs

Absolute value of first operand, so “4 abs” is equivalent to \(|4|\).

pow

First operand to the power of the second, so “4.3 5f pow” is equivalent to \(4.3^{5}\). Currently pow will only work on single or double precision floating point numbers or images. To be sure that a number is read as a floating point (even if it doesn’t have any non-zero decimals) put an f after it.

sqrt

The square root of the first operand, so “5 sqrt” is equivalent to \(\sqrt{5}\). The output type is determined from the input, so the output of this example will be 2 (since 5 doesn’t have any non-zero decimal digits). If you want 2.23607, run 5f sqrt instead, the f will ensure that a number will be read as a floating point number, even if it doesn’t have decimal digits. If the input image has an integer type, you should explicitly convert the image to floating point, for example a.fits float sqrt, see the type conversion operators below.

log

Natural logarithm of first operand, so “4 log” is equivalent to \(\ln(4)\). The output type is determined from the input, see the explanation under sqrt for more.

log10

Base-10 logarithm of first operand, so “4 log10” is equivalent to \(\log(4)\). The output type is determined from the input, see the explanation under sqrt for more.

minvalue

Minimum (non-blank) value in the top operand on the stack, so “a.fits minvalue” will push the the minimum pixel value in this image onto the stack. Therefore this operator is mainly intended for data (for example images), if the top operand is a number, this operator just returns it without any change. So note that when this operator acts on a single image, the output will no longer be an image, but a number. The output of this operand is in the same type as the input.

maxvalue

Maximum (non-blank) value of first operand in the same type, similar to minvalue.

numvalue

Number of non-blank elements in first operand in the uint64 type, similar to minvalue.

sumvalue

Sum of non-blank elements in first operand in the float32 type, similar to minvalue.

meanvalue

Mean value of non-blank elements in first operand in the float32 type, similar to minvalue.

stdvalue

Standard deviation of non-blank elements in first operand in the float32 type, similar to minvalue.

medianvalue

Median of non-blank elements in first operand with the same type, similar to minvalue.

min

The first popped operand to this operator must be a positive integer number which specifies how many further operands should be popped from the stack. The given number of operands must have the same type and size. Each pixel of the output of this operator will be set to the minimum value of the given number of operands (images) in that pixel.

For example the following command will produce an image with the same size and type as the inputs but each output pixel is set to the minimum respective pixel value of the three input images.

$ astarithmetic a.fits b.fits c.fits 3 min

Important notes:

• NaN/blank pixels will be ignored, see Blank pixels.
• The output will have the same type as the inputs. This is natural for the min and max operators, but for other similar operators (for example sum, or average) the per-pixel operations will be done in double precision floating point and then stored back in the input type. Therefore, if the input was an integer, C’s internal type conversion will be used.

max

Similar to min, but the pixels of the output will contain the maximum of the respective pixels in all operands in the stack.

num

Similar to min, but the pixels of the output will contain the number of the respective non-blank pixels in all input operands.

sum

Similar to min, but the pixels of the output will contain the sum of the respective pixels in all input operands.

mean

Similar to min, but the pixels of the output will contain the mean (average) of the respective pixels in all operands in the stack.

std

Similar to min, but the pixels of the output will contain the standard deviation of the respective pixels in all operands in the stack.

median

Similar to min, but the pixels of the output will contain the median of the respective pixels in all operands in the stack.

filter-mean

Apply mean filtering (or moving average) on the input dataset. During mean filtering, each pixel (data element) is replaced by the mean value of all its surrounding pixels (excluding blank values). The number of surrounding pixels in each dimension (to calculate the mean) is determined through the earlier operands that have been pushed onto the stack prior to the input dataset. The number of necessary operands is determined by the dimensions of the input dataset (first popped operand). The order of the dimensions on the command-line is the order in FITS format. Here is one example:

$ astarithmetic 5 4 image.fits filter-mean

In this example, each pixel is replaced by the mean of a 5 by 4 box around it. The box is 5 pixels along the first FITS dimension (horizontal when viewed in ds9) and 4 pixels along the second FITS dimension (vertical).

Each pixel will be placed in the center of the box that the mean is calculated on. If the given width along a dimension is even, then the center is assumed to be between the pixels (not in the center of a pixel). When the pixel is close to the edge, the pixels of the box that fall outside the image are ignored. Therefore, on the edge, less points will be used in calculating the mean.

The final effect of mean filtering is to smooth the input image, it is essentially a convolution with a kernel that has identical values for all its pixels (is flat), see Convolution process.

Note that blank pixels will also be affected by this operator: if there are any non-blank elements in the box surrounding a blank pixel, in the filtered image, it will have the mean of the non-blank elements, therefore it won’t be blank any more. If blank elements are important for your analysis, you can use the isblank with the where operator to set them back to blank after filtering.

filter-median

Apply median filtering on the input dataset. This is very similar to filter-mean, except that instead of the mean value of the box pixels, the median value is used to replace a pixel value. For more on how to use this operator, please see filter-mean.

The median is less susceptible to outliers compared to the mean. As a result, after median filtering, the pixel values will be more discontinuous than mean filtering.

filter-sigclip-mean

Apply a \(\sigma\)-clipped mean filtering onto the input dataset. This is very similar to filter-mean, except that all outliers (identified by the \(\sigma\)-clipping algorithm) have been removed, see Sigma clipping for more on the basics of this algorithm. As described there, two extra input parameters are necessary for \(\sigma\)-clipping: the multiple of \(\sigma\) and the termination criteria. filter-sigclip-mean therefore needs to pop two other operands from the stack after the dimensions of the box.

For example the line below uses the same box size as the example of filter-mean. However, all elements in the box that are iteratively beyond \(3\sigma\) of the distribution’s median are removed from the final calculation of the mean until the change in \(\sigma\) is less than \(0.2\).

$ astarithmetic 3 0.2 5 4 image.fits filter-sigclip-mean

The median (which needs a sorted dataset) is necessary for \(\sigma\)-clipping, therefore filter-sigclip-mean can be significantly slower than filter-mean. However, if there are strong outliers in the dataset that you want to ignore (for example emission lines on a spectrum when finding the continuum), this is a much better solution.

filter-sigclip-median

Apply a \(\sigma\)-clipped median filtering onto the input dataset. This operator and its necessary operands are almost identical to filter-sigclip-mean, except that after \(\sigma\)-clipping, the median value (which is less affected by outliers than the mean) is added back to the stack.

interpolate-medianngb

Interpolate all the blank elements of the second popped operand with the median of its nearest non-blank neighbors. The number of the nearest non-blank neighbors used to calculate the median is given by the first popped operand. Note that the distance of the nearest non-blank neighbors is irrelevant in this interpolation.

collapse-sum

Collapse the given dataset (second popped operand), by summing all elements along the first popped operand (a dimension in FITS standard: counting from one, from fastest dimension). The returned dataset has one dimension less compared to the input.

The output will have a double-precision floating point type irrespective of the input dataset’s type. Doing the operation in double-precision (64-bit) floating point will help the collapse (summation) be affected less by floating point errors. But afterwards, single-precision floating points are usually enough in real (noisy) datasets. So depending on the type of the input and its nature, it is recommended to use one of the type conversion operators on the returned dataset.

If any WCS is present, the returned dataset will also lack the respective dimension in its WCS matrix. Therefore, when the WCS is important for later processing, be sure that the input is aligned with the respective axises: all non-diagonal elements in the WCS matrix are zero.

One common application of this operator is the creation of pseudo broad-band or narrow-band 2D images from 3D data cubes. For example integral field unit (IFU) data products that have two spatial dimensions (first two FITS dimensions) and one spectral dimension (third FITS dimension). The command below will collapse the whole third dimension into a 2D array the size of the first two dimensions, and then convert the output to single-precision floating point (as discussed above).

$ astarithmetic cube.fits 3 collapse-sum float32

collapse-mean

Similar to collapse-sum, but the returned dataset will be the mean value along the collapsed dimension, not the sum.

collapse-number

Similar to collapse-sum, but the returned dataset will be the number of non-blank values along the collapsed dimension. The output will have a 32-bit signed integer type. If the input dataset doesn’t have blank values, all the elements in the returned dataset will have a single value (the length of the collapsed dimension). Therefore this is mostly relevant when there are blank values in the dataset.

collapse-min

Similar to collapse-sum, but the returned dataset will have the same numeric type as the input and will contain the minimum value for each pixel along the collapsed dimension.

collapse-max

Similar to collapse-sum, but the returned dataset will have the same numeric type as the input and will contain the maximum value for each pixel along the collapsed dimension.

erode

Erode the foreground pixels (with value 1) of the input dataset (second popped operand). The first popped operand is the connectivity (see description in connected-components). Erosion is simply a flipping of all foreground pixels (to background; with value 0) that are “touching” background pixels. “Touching” is defined by the connectivity. In effect, this carves off the outer borders of the foreground, making them thinner. This operator assumes a binary dataset (all pixels are 0 and 1).

dilate

Dilate the foreground pixels (with value 1) of the input dataset (second popped operand). The first popped operand is the connectivity (see description in connected-components). Erosion is simply a flipping of all background pixels (with value 0) to foreground that are “touching” foreground pixels. “Touching” is defined by the connectivity. In effect, this expands the outer borders of the foreground. This operator assumes a binary dataset (all pixels are 0 and 1).

connected-components

Find the connected components in the input dataset (second popped operand). The first popped is the connectivity used in the connected components algorithm. The second popped operand is the dataset where connected components are to be found. It is assumed to be a binary image (with values of 0 or 1). It must have an 8-bit unsigned integer type which is the format produced by conditional operators. This operator will return a labeled dataset where the non-zero pixels in the input will be labeled with a counter (starting from 1).

The connectivity is a number between 1 and the number of dimensions in the dataset (inclusive). 1 corresponds to the weakest (symmetric) connectivity between elements and the number of dimensions the strongest. For example on a 2D image, a connectivity of 1 corresponds to 4-connected neighbors and 2 corresponds to 8-connected neighbors.

One example usage of this operator can be the identification of regions above a certain threshold, as in the command below. With this command, Arithmetic will first separate all pixels greater than 100 into a binary image (where pixels with a value of 1 are above that value). Afterwards, it will label all those that are connected.

$ astarithmetic in.fits 100 gt 2 connected-components

If your input dataset doesn’t have a binary type, but you know all its values are 0 or 1, you can use the uint8 operator (below) to convert it to binary.

fill-holes

Flip background (0) pixels surrounded by foreground (1) in a binary dataset. This operator takes two operands (similar to connected-components): the first popped operand is the connectivity (to define a hole) and the second is the binary (0 or 1 valued) dataset to fill holes in.

invert

Invert an unsigned integer dataset. This is the only operator that ignores blank values (which are set to be the maximum values in the unsigned integer types).

This is useful in cases where the target(s) has(have) been imaged in absorption as raw formats (which are unsigned integer types). With this option, the maximum value for the given type will be subtracted from each pixel value, thus “inverting” the image, so the target(s) can be treated as emission. This can be useful when the higher-level analysis methods/tools only work on emission (positive skew in the noise, not negative).

lt

Less than: If the second popped (or left operand in infix notation, see Reverse polish notation) value is smaller than the first popped operand, then this function will return a value of 1, otherwise it will return a value of 0. If both operands are images, then all the pixels will be compared with their counterparts in the other image. If only one operand is an image, then all the pixels will be compared with the the single value (number) of the other operand. Finally if both are numbers, then the output is also just one number (0 or 1). When the output is not a single number, it will be stored as an unsigned char type.

le

Less or equal: similar to lt (‘less than’ operator), but returning 1 when the second popped operand is smaller or equal to the first.

gt

Greater than: similar to lt (‘less than’ operator), but returning 1 when the second popped operand is greater than the first.

ge

Greater or equal: similar to lt (‘less than’ operator), but returning 1 when the second popped operand is larger or equal to the first.

eq

Equality: similar to lt (‘less than’ operator), but returning 1 when the two popped operands are equal (to double precision floating point accuracy).

ne

Non-Equality: similar to lt (‘less than’ operator), but returning 1 when the two popped operands are not equal (to double precision floating point accuracy).

and

Logical AND: returns 1 if both operands have a non-zero value and 0 if both are zero. Both operands have to be the same kind: either both images or both numbers.

or

Logical OR: returns 1 if either one of the operands is non-zero and 0 only when both operators are zero. Both operands have to be the same kind: either both images or both numbers.

not

Logical NOT: returns 1 when the operand is zero and 0 when the operand is non-zero. The operand can be an image or number, for an image, it is applied to each pixel separately.

isblank

Test for a blank value (see Blank pixels). In essence, this is very similar to the conditional operators: the output is either 1 or 0 (see the ‘less than’ operator above). The difference is that it only needs one operand. Because of the definition of a blank pixel, a blank value is not even equal to itself, so you cannot use the equal operator above to select blank pixels. See the “Blank pixels” box below for more on Blank pixels in Arithmetic.

where

Change the input (pixel) value where/if a certain condition holds. The conditional operators above can be used to define the condition. Three operands are required for where. The input format is demonstrated in this simplified example:

$ astarithmetic modify.fits binary.fits if-true.fits where

The value of any pixel in modify.fits that corresponds to a non-zero and non-blank pixel of binary.fits will be changed to the value of the same pixel in if-true.fits (this may also be a number). The 3rd and 2nd popped operands (modify.fits and binary.fits respectively, see Reverse polish notation) have to have the same dimensions/size. if-true.fits can be either a number, or have the same dimension/size as the other two.

The 2nd popped operand (binary.fits) has to have uint8 (or unsigned char in standard C) type (see Numeric data types). It is treated as a binary dataset (with only two values: zero and non-zero, hence the name binary.fits in this example). However, commonly you won’t be dealing with an actual FITS file of a condition/binary image. You will probably define the condition in the same run based on some other reference image and use the conditional and logical operators above to make a true/false (or one/zero) image for you internally. For example the case below:

$ astarithmetic in.fits reference.fits 100 gt new.fits where

In the example above, any of the in.fits pixels that has a value in reference.fits greater than 100, will be replaced with the corresponding pixel in new.fits. Effectively the reference.fits 100 gt part created the condition/binary image which was added to the stack (in memory) and later used by where. The command above is thus equivalent to these two commands:

$ astarithmetic reference.fits 100 gt --output=binary.fits
$ astarithmetic in.fits binary.fits new.fits where

Finally, the input operands are read and used independently, so you can use the same file more than once as any of the operands.

When the 1st popped operand to where (if-true.fits) is a single number, it may be a NaN value (or any blank value, depending on its type) like the example below (see Blank pixels). When the number is blank, it will be converted to the blank value of the type of the 3rd popped operand (in.fits). Hence, in the example below, all the pixels in reference.fits that have a value greater than 100, will become blank in the natural data type of in.fits (even though NaN values are only defined for floating point types).

$ astarithmetic in.fits reference.fits 100 gt nan where

bitand

Bitwise AND operator: only bits with values of 1 in both popped operands will get the value of 1, the rest will be set to 0. For example (assuming numbers can be written as bit strings on the command-line): 00101000 00100010 bitand will give 00100000. Note that the bitwise operators only work on integer type datasets.

bitor

Bitwise inclusive OR operator: The bits where at least one of the two popped operands has a 1 value get a value of 1, the others 0. For example (assuming numbers can be written as bit strings on the command-line): 00101000 00100010 bitand will give 00101010. Note that the bitwise operators only work on integer type datasets.

bitxor

Bitwise exclusive OR operator: A bit will be 1 if it differs between the two popped operands. For example (assuming numbers can be written as bit strings on the command-line): 00101000 00100010 bitand will give 00001010. Note that the bitwise operators only work on integer type datasets.

lshift

Bitwise left shift operator: shift all the bits of the first operand to the left by a number of times given by the second operand. For example (assuming numbers can be written as bit strings on the command-line): 00101000 2 lshift will give 10100000. This is equivalent to multiplication by 4. Note that the bitwise operators only work on integer type datasets.

rshift

Bitwise right shift operator: shift all the bits of the first operand to the right by a number of times given by the second operand. For example (assuming numbers can be written as bit strings on the command-line): 00101000 2 rshift will give 00001010. Note that the bitwise operators only work on integer type datasets.

bitnot

Bitwise not (more formally known as one’s complement) operator: flip all the bits of the popped operand (note that this is the only unary, or single operand, bitwise operator). In other words, any bit with a value of 0 is changed to 1 and vice-versa. For example (assuming numbers can be written as bit strings on the command-line): 00101000 bitnot will give 11010111. Note that the bitwise operators only work on integer type datasets/numbers.

uint8

Convert the type of the popped operand to 8-bit unsigned integer type (see Numeric data types). The internal conversion of C will be used.

int8

Convert the type of the popped operand to 8-bit signed integer type (see Numeric data types). The internal conversion of C will be used.

uint16

Convert the type of the popped operand to 16-bit unsigned integer type (see Numeric data types). The internal conversion of C will be used.

int16

Convert the type of the popped operand to 16-bit signed integer (see Numeric data types). The internal conversion of C will be used.

uint32

Convert the type of the popped operand to 32-bit unsigned integer type (see Numeric data types). The internal conversion of C will be used.

int32

Convert the type of the popped operand to 32-bit signed integer type (see Numeric data types). The internal conversion of C will be used.

uint64

Convert the type of the popped operand to 64-bit unsigned integer (see Numeric data types). The internal conversion of C will be used.

float32

Convert the type of the popped operand to 32-bit (single precision) floating point (see Numeric data types). The internal conversion of C will be used.

float64

Convert the type of the popped operand to 64-bit (double precision) floating point (see Numeric data types). The internal conversion of C will be used.

set-AAA

Set the characters after the dash (AAA in the case shown here) as a name for the first popped operand on the stack. The named dataset will be freed from memory as soon as it is no longer needed, or if the name is reset to refer to another dataset later in the command. This operator thus enables re-usability of a dataset without having to re-read it from a file every time it is necessary during a process. When a dataset is necessary more than once, this operator can thus help simplify reading/writing on the command-line (thus avoiding potential bugs), while also speeding up the processing.

Like all operators, this operator pops the top operand off of the main processing stack, but unlike other operands, it won’t add anything back to the stack immediately. It will keep the popped dataset in memory through a separate list of named datasets (not on the main stack). That list will be used to add/copy any requested dataset to the main processing stack when the name is called.

The name to give the popped dataset is part of the operator’s name. For example the set-a operator of the command below, gives the name “a” to the contents of image.fits. This name is then used instead of the actual filename to multiply the dataset by two.

$ astarithmetic image.fits set-a a 2 x

The name can be any string, but avoid strings ending with standard filename suffixes (for example .fits)⁹².

One example of the usefulness of this operator is in the where operator. For example, let’s assume you want to mask all pixels larger than 5 in image.fits (extension number 1) with a NaN value. Without setting a name for the dataset, you have to read the file two times from memory in a command like this:

$ astarithmetic image.fits image.fits 5 gt nan where -g1

But with this operator you can simply give image.fits the name i and simplify the command above to the more readable one below (which greatly helps when the filename is long):

$ astarithmetic image.fits set-i   i i 5 gt nan where

$ astarithmetic input.fits nan eq --output=all-zeros.fits

Read in other formats.

JavaScript license information

GNU Astronomy Utilities 0.7 manual, August 2018.