[gnuastro-commits] master 3ba7c68b: Library (arithmetic.h): new operator to write a constant

[Mohammad Akhlaghi]

branch: master
commit 3ba7c68b1d1d1566100f30e2e64d85d7f094fac8
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>

    Library (arithmetic.h): new operator to write a constant
    
    Until now, there was no easy way to set all the elements within an input
    dataset to a constant value. We needed to use conditionals and 'where' for
    example.
    
    With this commit, a new 'constant' operator has been added to arithmetic
    for this purpose.
---
 NEWS                      |  5 ++++
 doc/gnuastro.texi         | 62 +++++++++++++++++++++++++++++++++++++++++++++++
 lib/arithmetic.c          | 49 ++++++++++++++++++++++++++++++++++---
 lib/gnuastro/arithmetic.h |  2 +-
 4 files changed, 114 insertions(+), 4 deletions(-)

diff --git a/NEWS b/NEWS
index db24a44e..6fc6f1f3 100644
--- a/NEWS
+++ b/NEWS
@@ -20,6 +20,11 @@ See the end of the file for license conditions.
      - counter: similar to 'index', but integers start with 1.
      - indexonly: similar to 'index', but pops the top stack dataset.
      - counteronly: similar to 'counter', but pops the top stack dataset.
+     - constant: replace all elements of input image or column with the
+       given constant value (for example to set all pixels in an image to
+       NaN or zero), or to give identifiers to each source catalog when
+       adding rows from different catalogs, see the documentation for
+       complete examples.
      - counts-to-nanomaggy: convert counts to Nanomaggy (with fixed
        zeropoint of 22.5, used as the pixel units of many surveys like
        SDSS). This was suggested by Giulia Golini.
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index ca630051..1656de7d 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -14829,6 +14829,14 @@ $ asttable a.fits --catrowfile=b.fits --catrowhdu=1 \
                   -cRA,DEC,MAGNITUDE --output=allrows.fits
 @end example
 
+@cartouche
+@cindex Provenance
+@noindent
+@strong{Provenance of each row:} When merging rows from separate catalogs, it 
is important to keep track of the source catalog of each row (its provenance).
+To do this, you can use @option{--catrowfile} in combination with the 
@code{constant} operator and @ref{Column arithmetic}.
+For a working example of this scenario, see the example within the 
documentation of the @code{constant} operator in @ref{Building new dataset and 
stack management}.
+@end cartouche
+
 @cartouche
 @noindent
 @strong{How to avoid repetition when adding rows:} this option will simply add 
the rows of multiple tables into one, it does not check their contents!
@@ -18609,6 +18617,56 @@ To further extend the example, you can use any of the 
noise-making operators to
 $ astarithmetic 100 200 2 makenew 5 mknoise-sigma
 @end example
 
+@item constant
+Return an operand that will have a constant value (first popped operand) in 
all its elements.
+The number of elements is read from the second popped operand.
+The second popped operand is only used for its number of elements, its numeric 
data type, or its values are fully ignored and it is later freed.
+
+@cindex Provenance
+Here is one useful scenario for this operator in tables: you want to merge the 
objects/rows of some catalogs together, but you first want to give each source 
catalog a label/counter that distinguishes between the source of each rows in 
the merged/final catalog (using @ref{Invoking asttable}).
+The steps below show the the usage of this.
+
+@example
+## Add label 1 to the RA, Dec, magnitude and magnitude error
+## rows of the first catalog.
+$ asttable cat-1.fits -cRA,DEC,MAG,MAG_ERR \
+       -c'arith $1 1 constant' --output=tab-1.fits
+
+## Similar to above, but for the second catalog, we don't
+## need metadata here (because '--catrowfile' takes the column
+## column meta data from its first input.
+$ asttable cat-2.fits -cRA,DEC,MAG,MAG_ERR \
+           -c'arith $1 2 constant' --output=tab-2.fits
+
+## Concatenate (merge/blend) the rows of the two tables into
+## one for the 5 columns, but also add a counter for each
+## object or row in the final catalog.
+$ asttable tab-1.fits --catrowfile=tab-2.fits \
+           -c'arith $1 counteronly' \
+           -cRA,DEC,MAG,MAG_ERR,5 --output=merged.fits \
+           --colmetadata=1,ID_MERGED,counter,"Merged ID." \
+           --colmetadata=6,SOURCE-CAT,counter,"Source ID."
+
+## Add keyword information on each input. It is very important
+## to preserve this within the merged catalog. If the tables
+## came from public databases (for example on VizieR), give
+## their public identifier as the value.
+$ astfits merged.fits --write=/,"Source catalogs" \
+          --write=CATSRC1,"I/355/gaiadr3","VizieR ID." \
+          --write=CATSRC2,"Jane Doe","Name of source."
+
+## Check the metadata in 'merged.fits' and clean the
+## temporary files.
+$ astfits merged.fits -h1
+@end example
+
+Like most operators, @code{constant} is not limited to tables, you can also 
apply it on images.
+In the example below, we'll use @code{constant} to set all the pixels of the 
input image to NaN (which is necessary in scenarios that you need to include in 
an image in an analysis, but you don't want its pixels to affect the 
processing):
+
+@example
+$ astarithmetic image.fits nan constant
+@end example
+
 @item swap
 Swap the top two operands on the stack.
 For example the @code{index} operator doesn't pop with the top operand (the 
input to @code{index}), it just adds the index image to the stack.
@@ -35274,6 +35332,10 @@ The number of dimensions is derived from the number of 
nodes in the list and the
 Just note that the list should be in the reverse of the desired dimensions.
 @end deffn
 
+@deffn  Macro GAL_ARITHMETIC_OP_MAKENEW
+Given a dataset and a constant,
+@end deffn
+
 @deffn Macro GAL_ARITHMETIC_OPSTR_LOADCOL_HDU
 @deffnx Macro GAL_ARITHMETIC_OPSTR_LOADCOL_FILE
 @deffnx Macro GAL_ARITHMETIC_OPSTR_LOADCOL_PREFIX
diff --git a/lib/arithmetic.c b/lib/arithmetic.c
index 92c77cb2..90ce54ba 100644
--- a/lib/arithmetic.c
+++ b/lib/arithmetic.c
@@ -2624,7 +2624,7 @@ arithmetic_box_around_ellipse(gal_data_t *d1, gal_data_t 
*d2,
 
 
 static gal_data_t *
-arithmetic_constant(int operator)
+arithmetic_constants_standard(int operator)
 {
   size_t one=1;
 
@@ -2814,6 +2814,41 @@ arithmetic_index_counter(gal_data_t *input, int 
operator, int flags)
 
 
 
+static gal_data_t *
+arithmetic_constant(gal_data_t *input, gal_data_t *constant, int operator,
+                    int flags)
+{
+  gal_data_t *out;
+  size_t i, tsize;
+
+  /* Sanity check. */
+  if(constant->size!=1)
+    error(EXIT_FAILURE, 0, "%s: the constant operand should only contain "
+          "a single element", __func__);
+
+  /* Allocate the necessary dataset. */
+  out = ( input->type==constant->type
+          ? input
+          : gal_data_alloc(NULL, constant->type, input->ndim, input->dsize,
+                           NULL, 0, input->minmapsize, input->quietmmap,
+                           NULL, NULL, NULL) );
+
+  /* Fill the output with the constant's value. */
+  tsize=gal_type_sizeof(out->type);
+  for(i=0;i<input->size;++i)
+    memcpy(gal_pointer_increment(out->array, i, out->type),
+           constant->array, tsize);
+
+  /* Clean up and return. */
+  if(out!=input && (flags & GAL_ARITHMETIC_FLAG_FREE))
+    gal_data_free(input);
+  return out;
+}
+
+
+
+
+
 gal_data_t *
 gal_arithmetic_load_col(char *str, int searchin, int ignorecase,
                         size_t minmapsize, int quietmmap)
@@ -3201,9 +3236,11 @@ gal_arithmetic_set_operator(char *string, size_t 
*num_operands)
   else if (!strcmp(string, "counteronly"))
     { op=GAL_ARITHMETIC_OP_COUNTERONLY;       *num_operands=1;  }
   else if (!strcmp(string, "makenew"))
-    { op=GAL_ARITHMETIC_OP_MAKENEW;           *num_operands=-1;  }
+    { op=GAL_ARITHMETIC_OP_MAKENEW;           *num_operands=-1; }
   else if (!strcmp(string, "size"))
     { op=GAL_ARITHMETIC_OP_SIZE;              *num_operands=2;  }
+  else if (!strcmp(string, "constant"))
+    { op=GAL_ARITHMETIC_OP_CONSTANT;          *num_operands=2;  }
 
   /* Operator not defined. */
   else
@@ -3346,6 +3383,7 @@ gal_arithmetic_operator_string(int operator)
 
     case GAL_ARITHMETIC_OP_SWAP:            return "swap";
     case GAL_ARITHMETIC_OP_INDEX:           return "index";
+    case GAL_ARITHMETIC_OP_CONSTANT:        return "constant";
     case GAL_ARITHMETIC_OP_COUNTER:         return "counter";
     case GAL_ARITHMETIC_OP_INDEXONLY:       return "indexonly";
     case GAL_ARITHMETIC_OP_COUNTERONLY:     return "counteronly";
@@ -3576,7 +3614,7 @@ gal_arithmetic(int operator, size_t numthreads, int 
flags, ...)
     case GAL_ARITHMETIC_OP_PI:
     case GAL_ARITHMETIC_OP_AVOGADRO:
     case GAL_ARITHMETIC_OP_FINESTRUCTURE:
-      out=arithmetic_constant(operator);
+      out=arithmetic_constants_standard(operator);
       break;
 
     /* Calculate the width and height of a box surrounding an ellipse with
@@ -3610,6 +3648,11 @@ gal_arithmetic(int operator, size_t numthreads, int 
flags, ...)
       d2 = va_arg(va, gal_data_t *);
       d1->next=d2; d2->next=NULL;  out=d1;
       break;
+    case GAL_ARITHMETIC_OP_CONSTANT:
+      d1 = va_arg(va, gal_data_t *);
+      d2 = va_arg(va, gal_data_t *);
+      out=arithmetic_constant(d1, d2, operator, flags);
+      break;
 
     /* When the operator is not recognized. */
     default:
diff --git a/lib/gnuastro/arithmetic.h b/lib/gnuastro/arithmetic.h
index 3e158519..7f5212c3 100644
--- a/lib/gnuastro/arithmetic.h
+++ b/lib/gnuastro/arithmetic.h
@@ -133,7 +133,6 @@ enum gal_arithmetic_operators
 
   GAL_ARITHMETIC_OP_E,            /* The base of natural logirithm. */
   GAL_ARITHMETIC_OP_PI,           /* Circle circumference by diameter. */
-
   GAL_ARITHMETIC_OP_C,            /* The speed of light. */
   GAL_ARITHMETIC_OP_G,            /* The gravitational constant. */
   GAL_ARITHMETIC_OP_H,            /* Plank's constant. */
@@ -211,6 +210,7 @@ enum gal_arithmetic_operators
 
   /* Meta operators */
   GAL_ARITHMETIC_OP_MAKENEW,      /* Build a new dataset, containing zeros.*/
+  GAL_ARITHMETIC_OP_CONSTANT,     /* Make a row with given constant. */
   GAL_ARITHMETIC_OP_SIZE,         /* Size of the dataset along an axis     */
   GAL_ARITHMETIC_OP_INDEX,        /* New with the index (counting from 0). */
   GAL_ARITHMETIC_OP_COUNTER,      /* New with the index (counting from 0). */