GNU Astronomy Utilities



6.2.4.19 Size and position operators

With the operators below you can get metadata about the top dataset on the stack.

index

Add a new operand to the stack with an integer type and the same size (in all dimensions) as top operand on the stack (before it was called; it is not popped!). The first pixel in the returned operand is zero, and every later pixel’s value is incremented by one. It is important to remember that the top operand is not popped by this operand, so it remains on the stack. After this operand is finished, it adds a new operand to the stack. To pop the previous operand, you can use the indexonly operator.

The data type of the output is always an unsigned integer, and its width is determined from the number of pixels/rows in the top operand. For example if there are only 108 rows in a table, the returned column will have an unsigned 8-bit integer type (that can keep 256 separate values). But if the top operand is a \(1000\times1000=10^6\) pixel image, the output will be a 32-bit unsigned integer. For the various types of integers, see Numeric data types.

To see the index image along with the actual image, you can use the --writeall operator to have a multi-HDU output (without --writeall, Arithmetic will complain if more than one operand is left at the end). After DS9 opens with the second command, flip between the two extensions.

$ astarithmetic image.fits index --writeall
$ astscript-fits-view image_arith.fits

Below is a review some usage examples of this operator:

Image: masking margins

With the command below, we will be masking all pixels that are 20 pixels away from the edges of the image (on the margin). Here is a description of the command below (for the basics of Arithmetic’s notation, see Reverse polish notation):

  • The index operator just adds a new dataset on the stack: unlike almost all other operators in Arithmetic, index doesn’t remove its input dataset from the stack (use indexonly for the “normal” behavior). This is because index returns the pixel metadata not data. As a result, after index, we have two operands on the stack: the input image and the index image.
  • With the set-i operator, the top operand (the image containing the index of each pixel) is popped from the stack and associated to the name i. Therefore after this, the stack only has the input image. For more on the set- operator, see Operand storage in memory or a file.
  • We need three values from the commands before Arithmetic (for the width and height of the image and the size of the margin). To make the rest of the command easier to read/use, we’ll define them in Arithmetic as three named operators (respectively called w, h and m). All three are integers that will have a positive value lower than \(2^{16}=65536\) (for a “normal” image!). Therefore, we will store them as 16-bit unsigned integers with the uint16 operator (this will help optimal processing in later steps). For more the type changing operators, see Numerical type conversion operators.
  • Using the modulo % and division (/) operators on the index image and the width, we extract the horizontal (X) and vertical (Y) positions of each pixel in separately named operands called X and Y. The maximum value in these two will also fit within an unsigned 16-bit integer, so we’ll also store these in that type.
  • For the horizontal (X) dimension, we select pixels that are less than the margin (X m lt) and those that are more than the width subtracted by the margin (X w m - gt).
  • The output of the lt and gt conditional operators above is a binary (0 or 1 valued) image. We therefore merge them into one binary image using the or operator. For more, see Conditional operators.
  • We repeat the two steps above for the vertical (Y) dimension.
  • Once the images containing the to-be-masked pixels in each dimension are made, we combine them into one binary image with a final or operator. At this point, the stack only has two operands: 1) the input image and 2) the binary image that has a value of 1 for all pixels whose value should be changed.
  • A single-element operand (nan) is added on the stack.
  • Using the where operator, we replace all the pixels that are non-zero in the second operand (on the margins) to the top operand’s value (NaN) in the third popped operand (image that was read from image.fits). For more on the where operator, see Conditional operators.
$ margin=20
$ width=$(astfits  image.fits --keyvalue=NAXIS1 -q)
$ height=$(astfits image.fits --keyvalue=NAXIS2 -q)
$ astarithmetic image.fits index       set-i \
                $width     uint16      set-w \
                $height    uint16      set-h \
                $margin    uint16      set-m \
                i w %      uint16      set-X \
                i w /      uint16      set-Y \
                X m lt     X w m - gt     or \
                Y m lt     Y h m - gt     or \
                or nan where
Image: Masking regions outside a circle

As another example for usage on an image, in the command below we are using index to define an image where each pixel contains the distance to the pixel with X,Y coordinates of 345,250. We are then using that distance image to only keep the pixels that are within a 50 pixel radius of that point.

The basic concept behind this process is very similar to the previous example, with a different mathematical definition for pixels to mask. The major difference is that we want the distance to a pixel within the image, we need to have negative values and the center coordinates can be in a sub-pixel positions. The best numeric datatype for intermediate steps is therefore floating point. 64-bit floating point can have a precision of up to 15 digits after the decimal point. This is far too much for what we need here: in astronomical imaging, the PSF is usually on the scale of 1 or more pixels (see Sampling theorem). So even reaching a precision of one millionth of a pixel (offered by 32-bit floating points) is beyond our wildest dreams (see Numeric data types). We will also define the horizontal (X) and vertical (Y) operands after shifting to the desired central point.

$ radius=50
$ centerx=345.2
$ centery=250.3
$ width=$(astfits image.fits --keyvalue=NAXIS1 -q)
$ astarithmetic image.fits index set-i \
                $width       uint16    set-w \
                $radius      float32   set-r \
                $centerx     float32   set-cx \
                $centery     float32   set-cy \
                i w % cx -             set-X \
                i w / cy -             set-Y \
                X X x Y Y x + sqrt r gt \
                nan where --output=arith-masked.fits

Optimal data types have significant benefits: choosing the minimum required datatype for your operation is very important to avoid wasting your CPU and RAM. Don’t simply default to 64-bit floating points for everything! Integer operations are much faster than floating points, and within floating point types, 32-bit is faster and will use half the RAM/storage! For more, see Numeric data types.

The example above was just a demo for usage of the index operator and some important concepts. But it is not the easiest way to achieve the desired result above! An easier way for the scenario above (to keep a circle within an image and set everything else to NaN) is to use MakeProfiles in combination with Arithmetic, like below:

$ radius=50
$ centerx=345.2
$ centery=250.3
$ echo "1 $centerx $centery 5 $radius 0 0 1 1 1" \
       | astmkprof --background=image.fits \
                   --mforflatpix --clearcanvas \
                   -omkprof-mask.fits --type=uint8
$ astarithmetic image.fits mkprof-mask.fits not \
                nan where -g1 -omkprof-masked.fits
Tables: adding new columns with row index

Within Table, you can use this operator to add an index column like below (see the counter operator for starting the count from one).

## The index will be the second column.
$ asttable table.fits -c'arith $1 index'

## The index will be the first column
$ asttable table.fits -c'arith $1 index swap'
indexonly

Similar to index, except that the top operand is popped from the stack and is no longer available afterwards.

counter

Similar to index, except that counting starts from one (not zero as in index). Counting from one is usually necessary when adding row counters in tables, like below:

$ asttable table.fits -c'arith $1 counter swap'
counteronly

Similar to counter, but the top operand before it is popped (no longer available).

size

Size of the dataset along a given FITS (or FORTRAN) dimension (counting from 1). The desired dimension should be the first popped operand and the dataset must be the second popped operand. The output will be a single unsigned integer (dimensions cannot be negative). For example, the following command will produce the size of the first extension/HDU (the default HDU) of a.fits along the second FITS axis.

$ astarithmetic a.fits 2 size

Not optimal: This operator reads the top element on the stack and then simply reads its size along the given dimension. On a small dataset this won’t consume much RAM, but if you want to put this in a pipeline or use it on large image, the extra RAM and slow operation can become meaningful. To avoid such issues, you can read the size along the given dimension using the --keyvalue option of Keyword inspection and manipulation. For example, in the code below, the X axis position of every pixel is returned:

$ width=$(astfits image.fits --keyvalue=NAXIS1 -q)
$ astarithmetic image.fits indexonly $width % -opix-x.fits