[gnuastro-commits] master d7a39359: Library (arithmetic): new zero point changing operator

[Mohammad Akhlaghi]

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

    Library (arithmetic): new zero point changing operator
    
    Until now, the only function that we had to convert zero points was to/from
    nanomaggy (zero point of 22.5). But this value is not suitable in some
    scenarios and other custom constant zero points may be necessary.
    
    With this commit, a new 'zeropoint-change' operator has been added to the
    arithmetic library (and thus the Arithmetic and Table programs) for this
    job. In order to do it, a new function in 'units.h' has been added to do
    the low-level operation.
---
 NEWS                      |  16 +++--
 doc/gnuastro.texi         |  27 ++++++++-
 lib/arithmetic.c          | 149 +++++++++++++++++++++++++++++++++++++++++-----
 lib/gnuastro/arithmetic.h |   1 +
 lib/gnuastro/units.h      |   4 ++
 lib/units.c               |  31 ++++++----
 6 files changed, 195 insertions(+), 33 deletions(-)

diff --git a/NEWS b/NEWS
index 1027cb77..c0e7da44 100644
--- a/NEWS
+++ b/NEWS
@@ -12,13 +12,18 @@ See the end of the file for license conditions.
     outputs as extra HDUs (by default, the output file is deleted). If the
     output does not exist, then this option has no affect.
 
-  - New operators:
+  - New operators (also avaialble in "Column Arithmetic" of the Table
+    progam):
 
     - filter-madclip-mean: filter/smooth the input using MAD-clipped mean.
 
     - filter-madclip-median: filter/smooth the input using MAD-clipped
       median.
 
+    - free: free (from memory) the top operand on the stack of
+      operands. This is useful in combination with operators that produce
+      more than one output operand.
+
     - madclip-maskfilled: mask (set to NaN) all input elements that are
       outliers (defined by MAD-clipping). Combined with the stacking
       operators this allows removing large contiugous outliers in your
@@ -27,9 +32,8 @@ See the end of the file for license conditions.
     - siglclip-maskfilled: similar to 'madclip-maskfilled', but defining
       outliers by Sigma-clipping.
 
-    - free: free (from memory) the top operand on the stack of
-      operands. This is useful in combination with operators that produce
-      more than one output operand.
+    - zeropoint-change: change the zero point of the input data set to a
+      new zero point.
 
 *** Statistics
 
@@ -88,6 +92,10 @@ See the end of the file for license conditions.
 *** Library
 - gal_statistics_concentration: measure the concentration of values around
   the median; see the book for the details.
+
+- gal_units_zeropoint_change: change the zero point of the input data set
+  to an output zero point.
+
 ** Removed features
 ** Changed features
 *** All programs
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 0e07ee13..f047122e 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -8935,8 +8935,8 @@ $ for f in i r g; do \
 @noindent
 @strong{Accounting for zero points:} An important step that we have not 
implemented in this section is to unify the zero points of the three filters.
 In the case of SDSS (and some other surveys), the images have already been 
brought to the same zero point, but that is not generally the case.
-So before subtracting sky (and estimating the standard deviation) you should 
also unify the zero points of your images (for example through Arithmetic's 
@code{counts-to-nanomaggy} or @code{counts-to-jy} described in @ref{Unit 
conversion operators}).
-If you don't already have the zero point of your images, see the dedicated 
tutorial: @ref{Zero point of an image}.
+So before subtracting sky (and estimating the standard deviation) you should 
also unify the zero points of your images (for example through Arithmetic's 
@code{counts-to-customzp}, @code{counts-to-nanomaggy} or @code{counts-to-jy} 
described in @ref{Unit conversion operators}).
+If you don't already have the zero point of your images, see the dedicated 
tutorial to measure it: @ref{Zero point of an image}.
 @end cartouche
 
 Now that we know the noise fluctuates around zero in all three images, we can 
start to define the values for the @option{--fluxlow} and @option{--fluxhigh}.
@@ -21867,6 +21867,18 @@ $ astarithmetic sdss-image.fits 22.5 counts-to-jy
 Convert Janskys to counts (usually CCD outputs) through an AB-magnitude based 
zero point.
 This is the inverse operation of the @code{counts-to-jy}, see there for usage 
example.
 
+@item zeropoint-change
+@cindex Zero point change
+Change the zero point of the input dataset to a new zero point.
+For example, with the command below, we are changing the zero point of an 
image from 26.892 to 27.0:
+
+@example
+$ astarithmetic image.fits 26.892 27.0 zeropoint-change
+@end example
+
+Note that the input or output zero points can also be an image (with the same 
size as the input image).
+This is very useful to correct situations where the zero point can change over 
the image (happens in single exposures) and you want to unify it across the 
whole image (to build a deep stacked image).
+
 @item counts-to-nanomaggy
 @cindex Nanomaggy
 Convert counts to Nanomaggy (with fixed zero point of 22.5, used as the pixel 
units of many surveys like SDSS).
@@ -21876,6 +21888,8 @@ For example if your image has a zero point of 24.93, 
you can convert it to Nanom
 $ astarithmetic image.fits 24.93 counts-to-nanomaggy
 @end example
 
+This is just a wrapper over the @code{zeropoint-change} operator where the 
output zero point is 22.5.
+
 @item nanomaggy-to-counts
 @cindex Nanomaggy
 Convert Nanomaggy to counts.
@@ -21886,6 +21900,8 @@ For example if you would like to convert an image in 
units of Nanomaggy (for exa
 $ astarithmetic image.fits 25.92 nanomaggy-to-counts
 @end example
 
+This is just a wrapper over the @code{zeropoint-change} operator where the 
input zero point is 22.5.
+
 @item mag-to-jy
 Convert AB magnitudes to Janskys, see @ref{Brightness flux magnitude}.
 
@@ -44517,14 +44533,21 @@ Convert the input value (assumed to be in 
Astronomical Units) to Parsecs.
 For the conversion equation, see the description of @code{au-to-pc} operator 
in @ref{Arithmetic operators}.
 @end deftypefun
 
+@deftypefun double gal_units_zeropoint_change (double @code{counts}, double 
@code{zeropoint_in}, double @code{zeropoint_out})
+@cindex Zero point change
+Convert the zero point of @code{counts} (which is @code{zeropoint_in}) to 
@code{zeropoint_out}.
+@end deftypefun
+
 @deftypefun double gal_units_counts_to_nanomaggy (double @code{counts}, double 
@code{zeropoint_ab})
 @cindex Nanomaggy
 @cindex Magnitude (nanomaggy)
 Convert counts to Nanomaggy (with fixed zero point of 22.5) through an AB 
magnitude-based zero point.
+This is just a wrapper around @code{gal_units_zeropoint_change} with 
@code{zeropoint_out=22.5}.
 @end deftypefun
 
 @deftypefun double gal_units_nanomaggy_to_counts (double @code{counts}, double 
@code{zeropoint_ab})
 Convert Nanomaggy (with fixed zero point of 22.5) to counts through an AB 
magnitude-based zero point.
+This is just a wrapper around @code{gal_units_zeropoint_change} with 
@code{zeropoint_in=22.5}.
 @end deftypefun
 
 @deftypefun double gal_units_pc_to_au (double @code{pc})
diff --git a/lib/arithmetic.c b/lib/arithmetic.c
index a073d54c..f5509f31 100644
--- a/lib/arithmetic.c
+++ b/lib/arithmetic.c
@@ -118,11 +118,36 @@ arithmetic_check_float_input(gal_data_t *in, int 
operator, char *numstr)
 */
 
 
+/* Given two datasets, if either one is a number, or both are the same size
+   (a good scenario) don't do anything. Otherwise, print an error
+   message. */
+static void
+arithmetic_sizes_bad_one_good(gal_data_t *a, gal_data_t *b, int flags,
+                              int operator)
+{
+  if( !( (flags & GAL_ARITHMETIC_FLAG_NUMOK) && (a->size==1 || b->size==1))
+      && gal_dimension_is_different(a, b) )
+    error(EXIT_FAILURE, 0, "%s: the non-number inputs to '%s' don't "
+          "have the same dimension/size", __func__,
+          gal_arithmetic_operator_string(operator));
+  return;
+}
 
 
 
 
 
+static gal_data_t *
+arithmetic_to_certain_type(gal_data_t *in, uint8_t desired_type, int flags)
+{
+  return ( in->type==desired_type
+           ? in
+           : gal_data_copy_to_new_type(in, desired_type) );
+}
+
+
+
+
 
 
 
@@ -2835,17 +2860,10 @@ arithmetic_binary(int operator, int flags, gal_data_t 
*l, gal_data_t *r)
     }
 
 
-  /* Simple sanity check on the input sizes. */
-  if( !( (flags & GAL_ARITHMETIC_FLAG_NUMOK) && (l->size==1 || r->size==1))
-      && gal_dimension_is_different(l, r) )
-    error(EXIT_FAILURE, 0, "%s: the non-number inputs to '%s' don't "
-          "have the same dimension/size", __func__,
-          gal_arithmetic_operator_string(operator));
-
-
-  /* Print a warning if the inputs are both integers, but have different
-     signs (the user needs to know that the output may not be what they
-     expect!).*/
+  /* Sanity checks: on the sizes and if the inputs are both integers, but
+     have different signs (the user needs to know that the output may not
+     be what they expect!).*/
+  arithmetic_sizes_bad_one_good(l, r, flags, operator);
   if( (flags & GAL_ARITHMETIC_FLAG_QUIET)==0 )
     arithmetic_binary_int_sanity_check(l, r, operator);
 
@@ -3012,11 +3030,8 @@ arithmetic_function_binary_flt(int operator, int flags, 
gal_data_t *il,
     }
 
 
-  /* Simple sanity check on the input sizes. */
-  if( !( (flags & GAL_ARITHMETIC_FLAG_NUMOK) && (il->size==1 || ir->size==1))
-      && gal_dimension_is_different(il, ir) )
-    error(EXIT_FAILURE, 0, "%s: the input datasets don't have the same "
-          "dimension/size", __func__);
+  /* Check on the input sizes (one can be a number). */
+  arithmetic_sizes_bad_one_good(il, ir, flags, operator);
 
 
   /* Convert the values to double precision floating point if they are
@@ -3128,6 +3143,98 @@ arithmetic_function_binary_flt(int operator, int flags, 
gal_data_t *il,
 
 
 
+/* Functions that take three arguments. */
+#define TERFUNC_RUN_FUNCTION(OP){                                       \
+    double *oa=o->array, *of=oa + o->size;                              \
+    double *la=l->array, *ma=m->array, *ra=r->array;                    \
+    do                                                                  \
+      {                                                                 \
+        *oa=OP(*la, *ma, *ra);                                          \
+        if(l->size>1) {++la;} if(m->size>1) {++ma;} if(r->size>1) {++ra;} \
+      }                                                                 \
+    while(++oa<of);                                                     \
+  }
+
+static gal_data_t *
+arithmetic_function_ternary_flt(int operator, int flags, gal_data_t *il,
+                               gal_data_t *im, gal_data_t *ir)
+
+{
+  struct wcsprm *owcs=NULL;
+  gal_data_t *l, *m, *r, *o=NULL;
+  size_t zero=0, ondim, *odsize, minmapsize;
+  int quietmmap=il->quietmmap && im->quietmmap && ir->quietmmap;
+
+  /* The datasets may be empty. In this case, the output should also be
+     empty (we can have tables and images with 0 rows or pixels!). */
+  if(    il->size==0 || il->array==NULL
+      || im->size==0 || im->array==NULL
+      || ir->size==0 || ir->array==NULL )
+    {
+      /* If we should free the inputs, then immediately free the middle
+         dataset which we do not need. Then check to see which one of the
+         inputs is already allocated and free it, then  */
+      if(flags & GAL_ARITHMETIC_FLAG_FREE) gal_data_free(il);
+      if(flags & GAL_ARITHMETIC_FLAG_FREE) gal_data_free(im);
+      if(flags & GAL_ARITHMETIC_FLAG_FREE) gal_data_free(ir);
+      o=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &zero, NULL, 0, -1, 1,
+                       NULL, NULL, NULL);
+    }
+
+  /* Simple sanity check on the input sizes (they should all either be a
+     single number, or have exactly the same number of dimensions). */
+  arithmetic_sizes_bad_one_good(il, im, flags, operator);
+  arithmetic_sizes_bad_one_good(il, ir, flags, operator);
+  arithmetic_sizes_bad_one_good(im, ir, flags, operator);
+
+  /* Convert the values to double precision floating point. */
+  l=arithmetic_to_certain_type(il, GAL_TYPE_FLOAT64, flags);
+  m=arithmetic_to_certain_type(im, GAL_TYPE_FLOAT64, flags);
+  r=arithmetic_to_certain_type(ir, GAL_TYPE_FLOAT64, flags);
+
+  /* Set the output parameters. */
+  minmapsize=(l->minmapsize < m->minmapsize
+              ? (l->minmapsize<r->minmapsize?l->minmapsize:r->minmapsize)
+              : (m->minmapsize<r->minmapsize?m->minmapsize:r->minmapsize));
+  if( l->size > m->size )
+    {
+      if( l->size > r->size ) { ondim=l->ndim; odsize=l->dsize; }
+      else                    { ondim=r->ndim; odsize=r->dsize; }
+    }
+  else
+    {
+      if( m->size > r->size ) { ondim=m->ndim; odsize=m->dsize; }
+      else                    { ondim=r->ndim; odsize=r->dsize; }
+    }
+  owcs = l->wcs ? l->wcs : (m->wcs ? m->wcs : r->wcs);
+
+  /* Allocate the output and do the operation */
+  o=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, ondim, odsize, owcs, 0,
+                   minmapsize, quietmmap, NULL, NULL, NULL);
+  switch(operator)
+    {
+    case GAL_ARITHMETIC_OP_ZEROPOINT_CHANGE:
+      TERFUNC_RUN_FUNCTION( gal_units_zeropoint_change ); break;
+    default:
+      error(EXIT_FAILURE, 0, "%s: operator code %d not recognized",
+            __func__, operator);
+    }
+
+  /* Clean up. */
+  if(l!=il) gal_data_free(l);
+  if(m!=im) gal_data_free(m);
+  if(r!=ir) gal_data_free(r);
+  if(flags & GAL_ARITHMETIC_FLAG_FREE)
+    { gal_data_free(il); gal_data_free(im); gal_data_free(ir); }
+
+  /* Return */
+  return o;
+}
+
+
+
+
+
 /* The list of arguments are:
      d1: Input data (counts or SB).
      d2: Zeropoint.
@@ -3979,6 +4086,8 @@ gal_arithmetic_set_operator(char *string, size_t 
*num_operands)
     { op=GAL_ARITHMETIC_OP_COUNTS_TO_JY;      *num_operands=2;  }
   else if (!strcmp(string, "jy-to-counts"))
     { op=GAL_ARITHMETIC_OP_JY_TO_COUNTS;      *num_operands=2;  }
+  else if (!strcmp(string, "zeropoint-change"))
+    { op=GAL_ARITHMETIC_OP_ZEROPOINT_CHANGE; *num_operands=3;}
   else if (!strcmp(string, "counts-to-nanomaggy"))
     { op=GAL_ARITHMETIC_OP_COUNTS_TO_NANOMAGGY; *num_operands=2;}
   else if (!strcmp(string, "nanomaggy-to-counts"))
@@ -4341,6 +4450,7 @@ gal_arithmetic_operator_string(int operator)
     case GAL_ARITHMETIC_OP_SB_TO_COUNTS:    return "sb-to-counts";
     case GAL_ARITHMETIC_OP_COUNTS_TO_JY:    return "counts-to-jy";
     case GAL_ARITHMETIC_OP_JY_TO_COUNTS:    return "jy-to-counts";
+    case GAL_ARITHMETIC_OP_ZEROPOINT_CHANGE: return "zeropoint-change";
     case GAL_ARITHMETIC_OP_COUNTS_TO_NANOMAGGY: return "counts-to-nanomaggy";
     case GAL_ARITHMETIC_OP_NANOMAGGY_TO_COUNTS: return "nanomaggy-to-counts";
     case GAL_ARITHMETIC_OP_MAG_TO_JY:       return "mag-to-jy";
@@ -4643,6 +4753,13 @@ gal_arithmetic(int operator, size_t numthreads, int 
flags, ...)
       out=arithmetic_function_binary_flt(operator, flags, d1, d2);
       break;
 
+    case GAL_ARITHMETIC_OP_ZEROPOINT_CHANGE:
+      d1 = va_arg(va, gal_data_t *);
+      d2 = va_arg(va, gal_data_t *);
+      d3 = va_arg(va, gal_data_t *);
+      out=arithmetic_function_ternary_flt(operator, flags, d1, d2, d3);
+      break;
+
     /* More complex operators. */
     case GAL_ARITHMETIC_OP_COUNTS_TO_SB:
     case GAL_ARITHMETIC_OP_SB_TO_COUNTS:
diff --git a/lib/gnuastro/arithmetic.h b/lib/gnuastro/arithmetic.h
index c5cc75b5..c8fb1ff2 100644
--- a/lib/gnuastro/arithmetic.h
+++ b/lib/gnuastro/arithmetic.h
@@ -157,6 +157,7 @@ enum gal_arithmetic_operators
   GAL_ARITHMETIC_OP_JY_TO_COUNTS, /* Janskys to Counts (AB zeropoint). */
   GAL_ARITHMETIC_OP_MAG_TO_JY,    /* Magnitude to Janskys. */
   GAL_ARITHMETIC_OP_JY_TO_MAG,    /* Janskys to Magnitude. */
+  GAL_ARITHMETIC_OP_ZEROPOINT_CHANGE, /* Change the zero point. */
   GAL_ARITHMETIC_OP_COUNTS_TO_NANOMAGGY,/* Counts to SDSS nanomaggies. */
   GAL_ARITHMETIC_OP_NANOMAGGY_TO_COUNTS,/* SDSS nanomaggies to counts. */
   GAL_ARITHMETIC_OP_AU_TO_PC,     /* Astronomical units (AU) to Parsecs (PC). 
*/
diff --git a/lib/gnuastro/units.h b/lib/gnuastro/units.h
index 4c06c082..c72d4c13 100644
--- a/lib/gnuastro/units.h
+++ b/lib/gnuastro/units.h
@@ -99,6 +99,10 @@ gal_units_counts_to_jy(double counts, double zeropoint_ab);
 double
 gal_units_jy_to_counts(double jy, double zeropoint_ab);
 
+double
+gal_units_zeropoint_change(double counts, double zeropoint_ab,
+                           double custom_zp);
+
 double
 gal_units_counts_to_nanomaggy(double counts, double zeropoint_ab);
 
diff --git a/lib/units.c b/lib/units.c
index e727f00a..2e255852 100644
--- a/lib/units.c
+++ b/lib/units.c
@@ -463,19 +463,30 @@ gal_units_jy_to_counts(double jy, double zeropoint_ab)
 
 
 
-/* Convert counts to nanomaggy. The job of this function is equivalent to
-   the double-call bellow. We just don't want to repeat some extra
-   multiplication/divisions.
+/* Convert counts to a custom zero point. The job of this function is
+   equivalent to the double-call bellow. We just don't want to repeat some
+   extra multiplication/divisions.
 
-     gal_units_jy_to_counts(gal_units_counts_to_jy(counts, zeropoint_ab),
-                            22.5)
+     gal_units_jy_to_counts(gal_units_counts_to_jy(counts, zeropoint_in),
+                            custom_out)
 */
 double
-gal_units_counts_to_nanomaggy(double counts, double zeropoint_ab)
+gal_units_zeropoint_change(double counts, double zeropoint_in,
+                           double zeropoint_out)
 {
   return ( counts
-           * pow(10, -1 * zeropoint_ab / 2.5)
-           / pow(10, -1 * 22.5         / 2.5) );
+           * pow(10, -1 * zeropoint_in  / 2.5)
+           / pow(10, -1 * zeropoint_out / 2.5) );
+}
+
+
+
+
+
+double
+gal_units_counts_to_nanomaggy(double counts, double zeropoint_ab)
+{
+  return gal_units_zeropoint_change(counts, zeropoint_ab, 22.5);
 }
 
 
@@ -485,9 +496,7 @@ gal_units_counts_to_nanomaggy(double counts, double 
zeropoint_ab)
 double
 gal_units_nanomaggy_to_counts(double counts, double zeropoint_ab)
 {
-  return ( counts
-           / pow(10, -1 * zeropoint_ab / 2.5)
-           * pow(10, -1 * 22.5         / 2.5) );
+  return gal_units_zeropoint_change(counts, 22.5, zeropoint_ab);
 }