GNU Astronomy Utilities


Next: , Previous: , Up: Gnuastro library   [Contents][Index]


12.3.14 Arithmetic on datasets (arithmetic.h)

When the dataset’s type and other information are already known, any programming language (including C) provides some very good tools for various operations (including arithmetic operations like addition) on the dataset with a simple loop. However, as an author of a program, making assumptions about the type of data, its dimensions and other basic characteristics will come with a large processing burden.

for example, if you always read your data as double precision floating points for a simple operation like addition with an integer constant, you will be wasting a lot of CPU and memory when the input dataset is int32 type for example, (see Numeric data types). This overhead may be small for small images, but as you scale your process up and work with hundred/thousands of files that can be very large, this overhead will take a significant portion of the processing power. The functions and macros in this section are designed precisely for this purpose: to allow you to do any of the defined operations on any dataset with no overhead (in the native type of the dataset).

Gnuastro’s Arithmetic program uses the functions and macros of this section, so please also have a look at the Arithmetic program and in particular Arithmetic operators for a better description of the operators discussed here.

The main function of this library is gal_arithmetic that is described below. It can take an arbitrary number of arguments as operands (depending on the operator, similar to printf). Its first two arguments are integers specifying the flags and operator. So first we will review the constants for the recognized flags and operators and discuss them, then introduce the actual function.

Macro: GAL_ARITHMETIC_FLAG_INPLACE
Macro: GAL_ARITHMETIC_FLAG_FREE
Macro: GAL_ARITHMETIC_FLAG_NUMOK
Macro: GAL_ARITHMETIC_FLAG_ENVSEED
Macro: GAL_ARITHMETIC_FLAG_QUIET
Macro: GAL_ARITHMETIC_FLAGS_BASIC

Bit-wise flags to pass onto gal_arithmetic (see below). To pass multiple flags, use the bitwise-or operator. for example, if you pass GAL_ARITHMETIC_FLAG_INPLACE | GAL_ARITHMETIC_FLAG_NUMOK, then the operation will be done in-place (without allocating a new array), and a single number will also be acceptable (that will be applied to all the pixels). Each flag is described below:

GAL_ARITHMETIC_FLAG_INPLACE

Do the operation in-place (in the input dataset, thus modifying it) to improve CPU and memory usage. If this flag is used, after gal_arithmetic finishes, the input dataset will be modified. It is thus useful if you have no more need for the input after the operation.

GAL_ARITHMETIC_FLAG_FREE

Free (all the) input dataset(s) after the operation is done. Hence the inputs are no longer usable after gal_arithmetic.

GAL_ARITHMETIC_FLAG_NUMOK

It is acceptable to use a number and an array together. for example, if you want to add all the pixels in an image with a single number you can pass this flag to avoid having to allocate a constant array the size of the image (with all the pixels having the same number).

GAL_ARITHMETIC_FLAG_ENVSEED

Use the pre-defined environment variable for setting the random number generator seed when an operator needs it (for example, mknoise-sigma). For more on random number generation in Gnuastro see Generating random numbers.

GAL_ARITHMETIC_FLAG_QUIET

Do Not print any warnings or messages for operators that may benefit from it. for example, by default the mknoise-sigma operator prints the random number generator function and seed that it used (in case the user wants to reproduce this result later). By activating this bit flag to the call, that extra information is not printed on the command-line.

GAL_ARITHMETIC_FLAGS_BASIC

A wrapper for activating the three “basic” operations that are commonly necessary together: GAL_ARITHMETIC_FLAG_INPLACE, GAL_ARITHMETIC_FLAG_FREE and GAL_ARITHMETIC_FLAG_NUMOK.

Macro: GAL_ARITHMETIC_OP_PLUS
Macro: GAL_ARITHMETIC_OP_MINUS
Macro: GAL_ARITHMETIC_OP_MULTIPLY
Macro: GAL_ARITHMETIC_OP_DIVIDE
Macro: GAL_ARITHMETIC_OP_LT
Macro: GAL_ARITHMETIC_OP_LE
Macro: GAL_ARITHMETIC_OP_GT
Macro: GAL_ARITHMETIC_OP_GE
Macro: GAL_ARITHMETIC_OP_EQ
Macro: GAL_ARITHMETIC_OP_NE
Macro: GAL_ARITHMETIC_OP_AND
Macro: GAL_ARITHMETIC_OP_OR

Binary operators (requiring two operands) that accept datasets of any recognized type (see Numeric data types). When gal_arithmetic is called with any of these operators, it expects two datasets as arguments. For a full description of these operators with the same name, see Arithmetic operators. The first dataset/operand will be put on the left of the operator and the second will be put on the right. The output type of the first four is determined from the input types (largest type of the inputs). The rest (which are all conditional operators) will output a binary uint8_t (or unsigned char) dataset with values of either 0 (zero) or 1 (one).

Macro: GAL_ARITHMETIC_OP_NOT

The logical NOT operator. When gal_arithmetic is called with this operator, it only expects one operand (dataset), since this is a unary operator. The output is uint8_t (or unsigned char) dataset of the same size as the input. Any non-zero element in the input will be 0 (zero) in the output and any 0 (zero) will have a value of 1 (one).

Macro: GAL_ARITHMETIC_OP_ISBLANK

A unary operator with output that is 1 for any element in the input that is blank, and 0 for any non-blank element. When gal_arithmetic is called with this operator, it will only expect one input dataset. The output dataset will have uint8_t (or unsigned char) type.

gal_arithmetic with this operator is just a wrapper for the gal_blank_flag function of Library blank values (blank.h) and this operator is just included for completeness in arithmetic operations. So in your program, it might be easier to just call gal_blank_flag.

Macro: GAL_ARITHMETIC_OP_WHERE

The three-operand where operator thoroughly discussed in Arithmetic operators. When gal_arithmetic is called with this operator, it will only expect three input datasets: the first (which is the same as the returned dataset) is the array that will be modified. The second is the condition dataset (that must have a uint8_t or unsigned char type), and the third is the value to be used if condition is non-zero.

As a result, note that the order of operands when calling gal_arithmetic with GAL_ARITHMETIC_OP_WHERE is the opposite of running Gnuastro’s Arithmetic program with the where operator (see Arithmetic). This is because the latter uses the reverse-Polish notation which is not necessary when calling a function (see Reverse polish notation).

Macro: GAL_ARITHMETIC_OP_SQRT
Macro: GAL_ARITHMETIC_OP_LOG
Macro: GAL_ARITHMETIC_OP_LOG10

Unary operator functions for calculating the square root (\(\sqrt{i}\)), \(ln(i)\) and \(log(i)\) mathematical operators on each element of the input dataset. The returned dataset will have a floating point type, but its precision is determined from the input: if the input is a 64-bit floating point, the output will also be 64-bit. Otherwise, the returned dataset will be 32-bit floating point: you do not gain precision by using these operators, but you gain in operating speed if you use the sufficient precision. See Numeric data types for more on the precision of floating point numbers to help in selecting your required floating point precision.

If you want your output to be 64-bit floating point but your input is a different type, you can convert the input to a 64-bit floating point type with gal_data_copy_to_new_type or gal_data_copy_to_new_type_free(see Copying datasets). Alternatively, you can use the GAL_ARITHMETIC_OP_TO_FLOAT64 operators in the arithmetic library.

Macro: GAL_ARITHMETIC_OP_SIN
Macro: GAL_ARITHMETIC_OP_COS
Macro: GAL_ARITHMETIC_OP_TAN
Macro: GAL_ARITHMETIC_OP_ASIN
Macro: GAL_ARITHMETIC_OP_ACOS
Macro: GAL_ARITHMETIC_OP_ATAN
Macro: GAL_ARITHMETIC_OP_ATAN2

Trigonometric functions (and their inverse). All the angles, either inputs or outputs, are in units of degrees.

Macro: GAL_ARITHMETIC_OP_SINH
Macro: GAL_ARITHMETIC_OP_COSH
Macro: GAL_ARITHMETIC_OP_TANH
Macro: GAL_ARITHMETIC_OP_ASINH
Macro: GAL_ARITHMETIC_OP_ACOSH
Macro: GAL_ARITHMETIC_OP_ATANH

Hyperbolic functions (and their inverse).

Macro: GAL_ARITHMETIC_OP_RA_TO_DEGREE
Macro: GAL_ARITHMETIC_OP_DEC_TO_DEGREE
Macro: GAL_ARITHMETIC_OP_DEGREE_TO_RA
Macro: GAL_ARITHMETIC_OP_DEGREE_TO_DEC

Unary operators to convert between degrees (as a single floating point number) to the sexagesimal Right Ascension and Declination format (as strings, respectively in the format of _h_m_s and _d_m_s). The first two operators expect a string operand (in the sexagesimal formats mentioned above, but also in the _:_:_) and will return a double-precision floating point operand. The latter two are the opposite.

Macro: GAL_ARITHMETIC_OP_COUNTS_TO_MAG
Macro: GAL_ARITHMETIC_OP_MAG_TO_COUNTS
Macro: GAL_ARITHMETIC_OP_MAG_TO_SB
Macro: GAL_ARITHMETIC_OP_SB_TO_MAG
Macro: GAL_ARITHMETIC_OP_COUNTS_TO_JY
Macro: GAL_ARITHMETIC_OP_JY_TO_COUNTS
Macro: GAL_ARITHMETIC_OP_MAG_TO_JY
Macro: GAL_ARITHMETIC_OP_JY_TO_MAG

Binary operators for converting brightness and surface brightness units to and from each other. The first operand to all of them are the values in the respective input unit. The second popped operand is the zero point, except for the operators that involve surface brightness (those with SB). For the surface brightness related operators, the second popped operand is the area in units of arcsec\(^2\).

Macro: GAL_ARITHMETIC_OP_COUNTS_TO_SB
Macro: GAL_ARITHMETIC_OP_SB_TO_COUNTS

Operators for converting counts to surface brightness and vice-versa. These operators take three operands: 1) the input dataset in units of counts or surface brightness (depending on the operator), 2) the zero point, 3) the area in units of arcsec\(^2\).

Macro: GAL_ARITHMETIC_OP_AU_TO_PC
Macro: GAL_ARITHMETIC_OP_PC_TO_AU
Macro: GAL_ARITHMETIC_OP_LY_TO_PC
Macro: GAL_ARITHMETIC_OP_PC_TO_LY
Macro: GAL_ARITHMETIC_OP_LY_TO_AU
Macro: GAL_ARITHMETIC_OP_AU_TO_LY

Unary operators to convert various distance units to and from each other: Astronomical Units (AU), Parsecs (PC) and Light years (LY).

Macro: GAL_ARITHMETIC_OP_MINVAL
Macro: GAL_ARITHMETIC_OP_MAXVAL
Macro: GAL_ARITHMETIC_OP_NUMBERVAL
Macro: GAL_ARITHMETIC_OP_SUMVAL
Macro: GAL_ARITHMETIC_OP_MEANVAL
Macro: GAL_ARITHMETIC_OP_STDVAL
Macro: GAL_ARITHMETIC_OP_MEDIANVAL

Unary operand statistical operators that will return a single value for datasets of any size. These are just wrappers around similar functions in Statistical operations (statistics.h) and are included in gal_arithmetic only for completeness (to use easily in Arithmetic). In your programs, it will probably be easier if you use those gal_statistics_ functions directly.

Macro: GAL_ARITHMETIC_OP_UNIQUE
Macro: GAL_ARITHMETIC_OP_NOBLANK

Unary operands that will remove some elements from the input dataset. The first will return the unique elements, and the second will return the non-blank elements. Due to the removal of elements, the dimensionality of the output will be lost.

These are just wrappers over the gal_statistics_unique and gal_blank_remove. These are just wrappers around similar functions in Statistical operations (statistics.h) and are included in gal_arithmetic only for completeness (to use easily in Arithmetic). In your programs, it will probably be easier if you use those gal_statistics_ functions directly.

Macro: GAL_ARITHMETIC_OP_ABS

Unary operand absolute-value operator.

Macro: GAL_ARITHMETIC_OP_MIN
Macro: GAL_ARITHMETIC_OP_MAX
Macro: GAL_ARITHMETIC_OP_NUMBER
Macro: GAL_ARITHMETIC_OP_SUM
Macro: GAL_ARITHMETIC_OP_MEAN
Macro: GAL_ARITHMETIC_OP_STD
Macro: GAL_ARITHMETIC_OP_MEDIAN

Multi-operand statistical operations. When gal_arithmetic is called with any of these operators, it will expect only a single operand that will be interpreted as a list of datasets (see List of gal_data_t). These operators can work on multiple threads using the numthreads argument. See the discussion under the min operator in Arithmetic operators.

The output will be a single dataset with each of its elements replaced by the respective statistical operation on the whole list. The type of the output is determined from the operator (irrespective of the input type): for GAL_ARITHMETIC_OP_MIN and GAL_ARITHMETIC_OP_MAX, it will be the same type as the input, for GAL_ARITHMETIC_OP_NUMBER, the output will be GAL_TYPE_UINT32 and for the rest, it will be GAL_TYPE_FLOAT32.

Macro: GAL_ARITHMETIC_OP_QUANTILE

Similar to the operands above (including GAL_ARITHMETIC_MIN), except that when gal_arithmetic is called with these operators, it requires two arguments. The first is the list of datasets like before, and the second is the 1-element dataset with the quantile value. The output type is the same as the inputs.

Macro: GAL_ARITHMETIC_OP_SIGCLIP_STD
Macro: GAL_ARITHMETIC_OP_SIGCLIP_MEAN
Macro: GAL_ARITHMETIC_OP_SIGCLIP_MEDIAN
Macro: GAL_ARITHMETIC_OP_SIGCLIP_NUMBER

Similar to the operands above (including GAL_ARITHMETIC_MIN), except that when gal_arithmetic is called with these operators, it requires two arguments. The first is the list of datasets like before, and the second is the 2-element list of \(\sigma\)-clipping parameters. The first element in the parameters list is the multiple of sigma and the second is the termination criteria (see Sigma clipping). The output type of GAL_ARITHMETIC_OP_SIGCLIP_NUMBER will be GAL_TYPE_UINT32 and for the rest it will be GAL_TYPE_FLOAT32.

Macro: GAL_ARITHMETIC_OP_MKNOISE_SIGMA
Macro: GAL_ARITHMETIC_OP_MKNOISE_POISSON
Macro: GAL_ARITHMETIC_OP_MKNOISE_UNIFORM

Add noise to the input dataset. These operators take two arguments: the first is the input data set (can have any dimensionality or number of elements. The second argument is the noise specifier (a single element, of any type): for a fixed-sigma noise, it is the Gaussian standard deviation, for the Poisson noise, it is the background (see Photon counting noise) and for the uniform distribution it is the width of the interval around each element of the input dataset.

By default, a separate random number generator seed will be used on each separate run of these operators. Therefore two identical runs on the same input will produce different results. You can get reproducible results by setting the GAL_RNG_SEED environment variable and activating the GAL_ARITHMETIC_FLAG_ENVSEED flag. For more on random number generation in Gnuastro, see Generating random numbers.

By default these operators will print the random number generator function and seed (in case the user wants to reproduce the result later), but this can be disabled by activating the bit-flag GAL_ARITHMETIC_FLAG_QUIET described above.

Macro: GAL_ARITHMETIC_OP_RANDOM_FROM_HIST
Macro: GAL_ARITHMETIC_OP_RANDOM_FROM_HIST_RAW

Select random values from a custom distribution (defined by a histogram). For more, see the description of the respective operators in Generating random numbers.

Macro: GAL_ARITHMETIC_OP_SIZE

Size operator that will return a single value for datasets of any kind. When gal_arithmetic is called with this operator, it requires two arguments. The first is the dataset, and the second is a single integer value. The output type is a single integer.

Macro: GAL_ARITHMETIC_OP_STITCH

Stitch a list of input datasets along the requested dimension. See the description of the stitch operator in Arithmetic (Dimensionality changing operators).

Macro: GAL_ARITHMETIC_OP_POW

Binary operator to-power operator. When gal_arithmetic is called with any of these operators, it will expect two operands: raising the first by the second (returning a floating point, inputs can be integers).

Macro: GAL_ARITHMETIC_OP_BITAND
Macro: GAL_ARITHMETIC_OP_BITOR
Macro: GAL_ARITHMETIC_OP_BITXOR
Macro: GAL_ARITHMETIC_OP_BITLSH
Macro: GAL_ARITHMETIC_OP_BITRSH
Macro: GAL_ARITHMETIC_OP_MODULO

Binary integer-only operand operators. These operators are only defined on integer data types. When gal_arithmetic is called with any of these operators, it will expect two operands: the first is put on the left of the operator and the second on the right. The ones starting with BIT are the respective bit-wise operators in C and MODULO is the modulo/remainder operator. For a discussion on these operators, please see Arithmetic operators.

The output type is determined from the input types and C’s internal conversions: it is strongly recommended that both inputs have the same type (any integer type), otherwise the bit-wise behavior will be determined by your compiler.

Macro: GAL_ARITHMETIC_OP_BITNOT

The unary bit-wise NOT operator. When gal_arithmetic is called with any of these operators, it will expect one operand of an integer type and preform the bitwise-NOT operation on it. The output will have the same type as the input.

Macro: GAL_ARITHMETIC_OP_TO_UINT8
Macro: GAL_ARITHMETIC_OP_TO_INT8
Macro: GAL_ARITHMETIC_OP_TO_UINT16
Macro: GAL_ARITHMETIC_OP_TO_INT16
Macro: GAL_ARITHMETIC_OP_TO_UINT32
Macro: GAL_ARITHMETIC_OP_TO_INT32
Macro: GAL_ARITHMETIC_OP_TO_UINT64
Macro: GAL_ARITHMETIC_OP_TO_INT64
Macro: GAL_ARITHMETIC_OP_TO_FLOAT32
Macro: GAL_ARITHMETIC_OP_TO_FLOAT64

Unary type-conversion operators. When gal_arithmetic is called with any of these operators, it will expect one operand and convert it to the requested type. Note that with these operators, gal_arithmetic is just a wrapper over the gal_data_copy_to_new_type or gal_data_copy_to_new_type_free that are discussed in Copying datasets. It accepts these operators only for completeness and easy usage in Arithmetic. So in your programs, it might be preferable to directly use those functions.

Macro: GAL_ARITHMETIC_OP_E
Macro: GAL_ARITHMETIC_OP_C
Macro: GAL_ARITHMETIC_OP_G
Macro: GAL_ARITHMETIC_OP_H
Macro: GAL_ARITHMETIC_OP_AU
Macro: GAL_ARITHMETIC_OP_LY
Macro: GAL_ARITHMETIC_OP_PI
Macro: GAL_ARITHMETIC_OP_AVOGADRO
Macro: GAL_ARITHMETIC_OP_FINESTRUCTURE

Return the respective mathematical constant. For their description please see Constants. The constant values are taken from the GNU Scientific Library’s headers (defined in gsl/gsl_math.h).

Macro: GAL_ARITHMETIC_OP_BOX_AROUND_ELLIPSE

Return the width (along horizontal) and height (along vertical) of a box that encompasses an ellipse with the same center point. For more on the three input operands to this operator see the description of box-around-ellipse. This function returns two datasets as a gal_data_t linked list. The top element of the list is the height and its next element is the width.

Macro: GAL_ARITHMETIC_OP_MAKENEW

Create a new, zero-valued dataset with an unsigned 8-bit data type. The length along each dimension of the dataset should be given as a single list of gal_data_ts. The number of dimensions is derived from the number of nodes in the list and the length along each dimension is the single-valued element within that list. Just note that the list should be in the reverse of the desired dimensions.

Macro: GAL_ARITHMETIC_OPSTR_LOADCOL_HDU
Macro: GAL_ARITHMETIC_OPSTR_LOADCOL_FILE
Macro: GAL_ARITHMETIC_OPSTR_LOADCOL_PREFIX
Macro: GAL_ARITHMETIC_OPSTR_LOADCOL_HDU_LEN
Macro: GAL_ARITHMETIC_OPSTR_LOADCOL_FILE_LEN
Macro: GAL_ARITHMETIC_OPSTR_LOADCOL_PREFIX_LEN

Constant components of the load-col- operator (see Loading external columns). These are just fixed strings (and their lengths) that are placed in between the various components of that operator to allow choosing a certain column of a certain HDU of a certain file.

Function:
gal_data_t *
gal_arithmetic (int operator, size_t numthreads, int flags, ...)

Apply the requested arithmetic operator on the operand(s). The operator is identified through the macros above (that start with GAL_ARITHMETIC_OP_). The number of necessary operands (number of arguments to replace ‘...’ in the declaration of this function, above) depends on the operator and is described under each operator, above. Each operand has a type of ‘gal_data_t *’ (see last paragraph with example).

If the operator can work on multiple threads, the number of threads can be specified with numthreads. When the operator is single-threaded, numthreads will be ignored. Special conditions can also be specified with the flag operator (a bit-flag with bits described above, for example, GAL_ARITHMETIC_FLAG_INPLACE or GAL_ARITHMETIC_FLAG_FREE).

gal_arithmetic is a multi-argument function (like C’s printf). In other words, the number of necessary arguments is not fixed and depends on the value to operator. Below, you can see a minimal, fully working example, showing how different operators need different numbers of arguments.

#include <stdio.h>
#include <stdlib.h>
#include <gnuastro/fits.h>
#include <gnuastro/arithmetic.h>

int
main(void)
{
  /* Define the datasets and flag. */
  gal_data_t *in1, *in2, *out1, *out2;
  int flag=GAL_ARITHMETIC_FLAGS_BASIC;

  /* Read the input images. */
  in1=gal_fits_img_read("image1.fits", "1", -1, 1);
  in2=gal_fits_img_read("image2.fits", "1", -1, 1);

  /* Take the logarithm (base-e) of the first input. */
  out1=gal_arithmetic(GAL_ARITHMETIC_OP_LOG, 1, flag, in1);

  /* Add the second input with the logarithm of the first. */
  out2=gal_arithmetic(GAL_ARITHMETIC_OP_PLUS, 1, flag, in2, out1);

  /* Write the output into a file. */
  gal_fits_img_write(out2, "out.fits", NULL, NULL);

  /* Clean up. Due to the in-place flag (in
   * 'GAL_ARITHMETIC_FLAGS_BASIC'), 'out1' and 'out2' point to the
   * same array in memory and due to the freeing flag, any input
   * dataset(s) that were not returned have been freed internally
   * by 'gal_arithmetic'. Therefore it is only necessary to free
   * 'out2': all other allocated spaces have been freed internally.
   * before reaching this point. */
  gal_data_free(out2);

  /* Return control back to the OS (saying that we succeeded). */
  return EXIT_SUCCESS;
}

As you see above, you can feed the returned dataset from one call of gal_arithmetic to another call. The advantage of using gal_arithmetic (as opposed to manually writing a for or while loop and doing the operation with the + operator and log() function yourself), is that you do not have to worry about the type of the input data (for a list of acceptable data types in Gnuastro, see Library data types (type.h)). Arithmetic will automatically deal with the data types internally and choose the best output type depending on the operator.

Function:
int
gal_arithmetic_set_operator (char *string, size_t *num_operands)

Return the operator macro/code that corresponds to string. The number of operands that it needs are written into the space that *num_operands points to. If the string could not be interpreted as an operator, this function will return GAL_ARITHMETIC_OP_INVALID.

This function will check string with the fixed human-readable names (using strcmp) for the operators and return the two numbers. Note that string must only contain the single operator name and nothing else (not even any extra white space).

Function:
char *
gal_arithmetic_operator_string (int operator)

Return the human-readable standard string that corresponds to the given operator. for example, when the input is GAL_ARITHMETIC_OP_PLUS or GAL_ARITHMETIC_OP_MEAN, the strings + or mean will be returned.

Function:
gal_data_t *
gal_arithmetic_load_col (char *str, int searchin, int ignorecase, size_t minmapsize, int quietmmap)

Return the column that corresponds to the identifier in the input string (str). str is expected to be in the format of the load-col- operator (see Loading external columns). This function will extract the column identifier, the file name and the HDU (if necessary) from the string, read the requested column in memory and return it.

See Table input output (table.h) for the macros that can be given to searchin and ignorecase and Generic data container (gal_data_t) for the definitions of minmapsize and quietmmap.


Next: Tessellation library (tile.h), Previous: World Coordinate System (wcs.h), Up: Gnuastro library   [Contents][Index]