[gnuastro-commits] master 4ef33c9: Arithmetic: coadding operators return float32 and uint32 (for counting)

[Mohammad Akhlaghi]

branch: master
commit 4ef33c9d0b8f1114ffdc8963878826f69ad163d3
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>

    Arithmetic: coadding operators return float32 and uint32 (for counting)
    
    Until now, the returned dataset from the coadding operators had the same
    type as the input datasets. This would cause a lot of information (for
    example when the inputs are integers and we want to coadd by mean, which
    needs float), or have too much information (for example when the inputs are
    float and you want to coadd by number, which needs integer).
    
    To solve this problem, with this commit, the output type is determined by
    the operator, not the input data type.
    
    This was suggested by Raul Infante-Sainz.
---
 doc/gnuastro.texi | 96 +++++++++++++++++++++++++++++++++++--------------------
 lib/arithmetic.c  | 42 +++++++++++++++++++-----
 2 files changed, 96 insertions(+), 42 deletions(-)

diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index b41a0ca..c75c218 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -12668,6 +12668,9 @@ to @command{minvalue}.
 
 @cindex NaN
 @item min
+For each pixel, find the minimum value in all given datasets. The output
+will have the same type as the input.
+
 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. All the subsequently popped operands must have the same type and
@@ -12703,37 +12706,50 @@ conversion will be used.
 @end itemize
 
 @item max
-Similar to @command{min}, but the pixels of the output will contain
-the maximum of the respective pixels in all operands in the stack.
+For each pixel, find the maximum value in all given datasets. The output
+will have the same type as the input. This operator is called similar to
+the @command{min} operator, please see there for more.
 
 @item number
-Similar to @command{min}, but the pixels of the output will contain the
-number of the respective non-blank pixels in all input operands.
+For each pixel count the number of non-blank pixels in all given
+datasets. The output will be an unsigned 32-bit integer datatype (see
address@hidden data types}). This operator is called similar to the
address@hidden operator, please see there for more.
 
 @item sum
-Similar to @command{min}, but the pixels of the output will contain the sum
-of the respective pixels in all input operands.
+For each pixel, calculate the sum in all given datasets. The output will
+have the a single-precision (32-bit) floating point type. This operator is
+called similar to the @command{min} operator, please see there for more.
 
 @item mean
-Similar to @command{min}, but the pixels of the output will contain the
-mean (average) of the respective pixels in all operands in the stack.
+For each pixel, calculate the mean in all given datasets. The output will
+have the a single-precision (32-bit) floating point type. This operator is
+called similar to the @command{min} operator, please see there for more.
 
 @item std
-Similar to @command{min}, but the pixels of the output will contain the
-standard deviation of the respective pixels in all operands in the stack.
+For each pixel, find the standard deviation in all given datasets. The
+output will have the a single-precision (32-bit) floating point type. This
+operator is called similar to the @command{min} operator, please see there
+for more.
 
 @item median
-Similar to @command{min}, but the pixels of the output will contain
-the median of the respective pixels in all operands in the stack.
+For each pixel, find the median in all given datasets. The output will have
+the a single-precision (32-bit) floating point type. This operator is
+called similar to the @command{min} operator, please see there for more.
 
 @item sigclip-number
-Combine the specified number of inputs into a single output that contains
-the number of remaining elements after @mymath{\sigma}-clipping on each
-element/pixel (for more on @mymath{\sigma}-clipping, see @ref{Sigma
-clipping}). This operator is very similar to @command{min}, with the
-exception that it expects two operands (parameters for sigma-clipping)
-before the total number of inputs. The first popped operand is the
-termination criteria and the second is the multiple of @mymath{\sigma}.
+For each pixel, find the sigma-clipped number (after removing outliers) in
+all given datasets. The output will have the an unsigned 32-bit integer
+type (see @ref{Numeric data types}).
+
+This operator will combine the specified number of inputs into a single
+output that contains the number of remaining elements after
address@hidden on each element/pixel (for more on
address@hidden, see @ref{Sigma clipping}). This operator is very
+similar to @command{min}, with the exception that it expects two operands
+(parameters for sigma-clipping) before the total number of inputs. The
+first popped operand is the termination criteria and the second is the
+multiple of @mymath{\sigma}.
 
 For example in the command below, the first popped operand (@command{0.2})
 is the sigma clipping termination criteria. If the termination criteria is
@@ -12749,19 +12765,22 @@ astarithmetic a.fits b.fits c.fits 3 5 0.2 
sigclip-number
 @end example
 
 @item sigclip-median
-Combine the specified number of inputs into a single output that contains
-the @emph{median} after @mymath{\sigma}-clipping on each element/pixel. For
-more see @command{sigclip-number}.
+For each pixel, find the sigma-clipped median in all given datasets. The
+output will have the a single-precision (32-bit) floating point type. This
+operator is called similar to the @command{sigclip-number} operator, please
+see there for more.
 
 @item sigclip-mean
-Combine the specified number of inputs into a single output that contains
-the @emph{mean} after @mymath{\sigma}-clipping on each element/pixel. For
-more see @command{sigclip-number}.
+For each pixel, find the sigma-clipped mean in all given datasets. The
+output will have the a single-precision (32-bit) floating point type. This
+operator is called similar to the @command{sigclip-number} operator, please
+see there for more.
 
 @item sigclip-std
-Combine the specified number of inputs into a single output that contains
-the @emph{standard deviation} after @mymath{\sigma}-clipping on each
-element/pixel. For more see @command{sigclip-number}.
+For each pixel, find the sigma-clipped standard deviation in all given
+datasets. The output will have the a single-precision (32-bit) floating
+point type. This operator is called similar to the @command{sigclip-number}
+operator, please see there for more.
 
 @item filter-mean
 Apply mean filtering (or @url{https://en.wikipedia.org/wiki/Moving_average,
@@ -27708,11 +27727,18 @@ Unary operand absolute-value operator.
 @deffnx Macro GAL_ARITHMETIC_OP_MEDIAN
 Multi-operand statistical operations. When @code{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 @ref{List of gal_data_t}. The
-output will be a single dataset with each of its elements replaced by the
-respective statistical operation on the whole list. These operators can
-work on multiple threads using the @code{numthreads} argument. See the
-discussion under the @code{min} operator in @ref{Arithmetic operators}.
+be interpreted as a list of datasets (see @ref{List of gal_data_t}). These
+operators can work on multiple threads using the @code{numthreads}
+argument. See the discussion under the @code{min} operator in
address@hidden 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 @code{GAL_ARITHMETIC_OP_MIN} and @code{GAL_ARITHMETIC_OP_MAX}, it will
+be the same type as the input, for @code{GAL_ARITHMETIC_OP_NUMBER}, the
+output will be @code{GAL_TYPE_UINT32} and for the rest, it will be
address@hidden
 @end deffn
 
 @deffn  Macro GAL_ARITHMETIC_OP_SIGCLIP_STD
@@ -27724,7 +27750,9 @@ that when @code{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 @ref{Sigma clipping}).
+the termination criteria (see @ref{Sigma clipping}). The output type of
address@hidden will be @code{GAL_TYPE_UINT32} and
+for the rest it will be @code{GAL_TYPE_FLOAT32}.
 @end deffn
 
 @deffn Macro GAL_ARITHMETIC_OP_POW
diff --git a/lib/arithmetic.c b/lib/arithmetic.c
index dfdfabb..9250588 100644
--- a/lib/arithmetic.c
+++ b/lib/arithmetic.c
@@ -686,8 +686,8 @@ struct multioperandparams
 
 
 #define MULTIOPERAND_MIN(TYPE) {                                        \
-    TYPE t, max;                                                        \
     size_t n, j=0;                                                      \
+    TYPE t, max, *o=p->out->array;                                      \
     gal_type_max(p->list->type, &max);                                  \
                                                                         \
     /* Go over all the pixels assigned to this thread. */               \
@@ -716,8 +716,8 @@ struct multioperandparams
 
 
 #define MULTIOPERAND_MAX(TYPE) {                                        \
-    TYPE t, min;                                                        \
     size_t n, j=0;                                                      \
+    TYPE t, min, *o=p->out->array;                                      \
     gal_type_min(p->list->type, &min);                                  \
                                                                         \
     /* Go over all the pixels assigned to this thread. */               \
@@ -748,6 +748,7 @@ struct multioperandparams
 #define MULTIOPERAND_NUM {                                              \
     int use;                                                            \
     size_t n, j;                                                        \
+    uint32_t *o=p->out->array;                                          \
                                                                         \
     /* Go over all the pixels assigned to this thread. */               \
     for(tind=0; tprm->indexs[tind] != GAL_BLANK_SIZE_T; ++tind)         \
@@ -780,6 +781,7 @@ struct multioperandparams
     int use;                                                            \
     double sum;                                                         \
     size_t n, j;                                                        \
+    float *o=p->out->array;                                             \
                                                                         \
     /* Go over all the pixels assigned to this thread. */               \
     for(tind=0; tprm->indexs[tind] != GAL_BLANK_SIZE_T; ++tind)         \
@@ -813,6 +815,7 @@ struct multioperandparams
     int use;                                                            \
     double sum;                                                         \
     size_t n, j;                                                        \
+    float *o=p->out->array;                                             \
                                                                         \
     /* Go over all the pixels assigned to this thread. */               \
     for(tind=0; tprm->indexs[tind] != GAL_BLANK_SIZE_T; ++tind)         \
@@ -846,6 +849,7 @@ struct multioperandparams
     int use;                                                            \
     size_t n, j;                                                        \
     double sum, sum2;                                                   \
+    float *o=p->out->array;                                             \
                                                                         \
     /* Go over all the pixels assigned to this thread. */               \
     for(tind=0; tprm->indexs[tind] != GAL_BLANK_SIZE_T; ++tind)         \
@@ -883,6 +887,7 @@ struct multioperandparams
 #define MULTIOPERAND_MEDIAN(TYPE, QSORT_F) {                            \
     int use;                                                            \
     size_t n, j;                                                        \
+    float *o=p->out->array;                                             \
     TYPE *pixs=gal_pointer_allocate(p->list->type, p->dnum, 0,          \
                                     __func__, "pixs");                  \
                                                                         \
@@ -926,9 +931,10 @@ struct multioperandparams
 
 
 #define MULTIOPERAND_SIGCLIP(TYPE) {                                    \
-    float *sarr;                                                        \
     size_t n, j;                                                        \
     gal_data_t *sclip;                                                  \
+    uint32_t *N=p->out->array;                                          \
+    float *sarr, *o=p->out->array;                                      \
     TYPE *pixs=gal_pointer_allocate(p->list->type, p->dnum, 0,          \
                                     __func__, "pixs");                  \
     gal_data_t *cont=gal_data_alloc(pixs, p->list->type, 1, &p->dnum,   \
@@ -954,7 +960,7 @@ struct multioperandparams
               case GAL_ARITHMETIC_OP_SIGCLIP_STD:    o[j]=sarr[3]; break;\
               case GAL_ARITHMETIC_OP_SIGCLIP_MEAN:   o[j]=sarr[2]; break;\
               case GAL_ARITHMETIC_OP_SIGCLIP_MEDIAN: o[j]=sarr[1]; break;\
-              case GAL_ARITHMETIC_OP_SIGCLIP_NUMBER: o[j]=sarr[0]; break;\
+              case GAL_ARITHMETIC_OP_SIGCLIP_NUMBER: N[j]=sarr[0]; break;\
               default:                                                  \
                 error(EXIT_FAILURE, 0, "%s: a bug! the code %d is not " \
                       "valid for sigma-clipping results", __func__,     \
@@ -979,9 +985,9 @@ struct multioperandparams
 
 
 #define MULTIOPERAND_TYPE_SET(TYPE, QSORT_F) {                          \
+    TYPE b, **a;                                                        \
     gal_data_t *tmp;                                                    \
     size_t i=0, tind;                                                   \
-    TYPE b, **a, *o=p->out->array;                                      \
                                                                         \
     /* Allocate space to keep the pointers to the arrays of each. */    \
     /* Input data structure. The operators will increment these */      \
@@ -1110,9 +1116,9 @@ static gal_data_t *
 arithmetic_multioperand(int operator, int flags, gal_data_t *list,
                         gal_data_t *params, size_t numthreads)
 {
-  uint8_t *hasblank;
   size_t i=0, dnum=1;
   float p1=NAN, p2=NAN;
+  uint8_t *hasblank, otype;
   struct multioperandparams p;
   gal_data_t *out, *tmp, *ttmp;
 
@@ -1160,11 +1166,31 @@ arithmetic_multioperand(int operator, int flags, 
gal_data_t *list,
     }
 
 
+  /* Set the output dataset type */
+  switch(operator)
+    {
+    case GAL_ARITHMETIC_OP_MIN:            otype=list->type;       break;
+    case GAL_ARITHMETIC_OP_MAX:            otype=list->type;       break;
+    case GAL_ARITHMETIC_OP_NUMBER:         otype=GAL_TYPE_UINT32;  break;
+    case GAL_ARITHMETIC_OP_SUM:            otype=GAL_TYPE_FLOAT32; break;
+    case GAL_ARITHMETIC_OP_MEAN:           otype=GAL_TYPE_FLOAT32; break;
+    case GAL_ARITHMETIC_OP_STD:            otype=GAL_TYPE_FLOAT32; break;
+    case GAL_ARITHMETIC_OP_MEDIAN:         otype=GAL_TYPE_FLOAT32; break;
+    case GAL_ARITHMETIC_OP_SIGCLIP_STD:    otype=GAL_TYPE_FLOAT32; break;
+    case GAL_ARITHMETIC_OP_SIGCLIP_MEAN:   otype=GAL_TYPE_FLOAT32; break;
+    case GAL_ARITHMETIC_OP_SIGCLIP_MEDIAN: otype=GAL_TYPE_FLOAT32; break;
+    case GAL_ARITHMETIC_OP_SIGCLIP_NUMBER: otype=GAL_TYPE_UINT32;  break;
+    default:
+      error(EXIT_FAILURE, 0, "%s: operator code %d isn't recognized",
+            __func__, operator);
+    }
+
+
   /* Set the output data structure. */
-  if(flags & GAL_ARITHMETIC_INPLACE)
+  if( (flags & GAL_ARITHMETIC_INPLACE) && otype==list->type)
     out = list;                 /* The top element in the list. */
   else
-    out = gal_data_alloc(NULL, list->type, list->ndim, list->dsize,
+    out = gal_data_alloc(NULL, otype, list->ndim, list->dsize,
                          list->wcs, 0, list->minmapsize, NULL, NULL, NULL);