[gnuastro-commits] master c4ae8f4: Collapsing by extremum in library and Arithmetic program

[Mohammad Akhlaghi]

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

    Collapsing by extremum in library and Arithmetic program
    
    In some contexts, it is useful to collapse a given dataset according to the
    minimum/maximum along the given dimension. With this commit, a new
    `gal_dimension_collapse_minmax' was added to the library and `collapse-min'
    and `collapse-max' (which use that library function) were added to the
    Arithmetic program.
---
 NEWS                        |   3 ++
 bin/arithmetic/arithmetic.c |  14 ++++++
 bin/arithmetic/arithmetic.h |   2 +
 doc/gnuastro.texi           |  19 ++++++++
 lib/dimension.c             | 105 ++++++++++++++++++++++++++++++++++++++++----
 lib/gnuastro/dimension.h    |   3 +-
 6 files changed, 136 insertions(+), 10 deletions(-)

diff --git a/NEWS b/NEWS
index 388ec6d..4705876 100644
--- a/NEWS
+++ b/NEWS
@@ -14,6 +14,8 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
                possibly using it many times.
     - `fill-holes': Flip background (0) pixels surrounded by foreground (1).
     - `collapse-sum': collapse/remove a dimension by summing over it.
+    - `collapse-min': collapse/remove a dimension by using maximum value.
+    - `collapse-max': collapse/remove a dimension by using minimum value.
     - `collapse-mean': collapse/remove a dimension by averaging over it.
     - `collapse-number': Number of elements included in the collapse.
 
@@ -31,6 +33,7 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
     - gal_dimension_collapse_sum: collapse/remove a dimension by summing.
     - gal_dimension_collapse_mean: collapse/remove a dimension by averaging.
     - gal_dimension_collapse_number: collapse/remove a dimension by number.
+    - gal_dimension_collapse_minmax: collapse/remove a dimension by extremum.
     - gal_wcs_remove_dimension: Remove a dimension in the given WCS structure.
 
 ** Removed features
diff --git a/bin/arithmetic/arithmetic.c b/bin/arithmetic/arithmetic.c
index c2a2bb1..e511e69 100644
--- a/bin/arithmetic/arithmetic.c
+++ b/bin/arithmetic/arithmetic.c
@@ -751,6 +751,14 @@ arithmetic_collapse(struct arithmeticparams *p, char 
*token, int operator)
       collapsed=gal_dimension_collapse_number(input, input->ndim-dim);
       break;
 
+    case ARITHMETIC_OP_COLLAPSE_MIN:
+      collapsed=gal_dimension_collapse_minmax(input, input->ndim-dim, 0);
+      break;
+
+    case ARITHMETIC_OP_COLLAPSE_MAX:
+      collapsed=gal_dimension_collapse_minmax(input, input->ndim-dim, 1);
+      break;
+
     default:
       error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix the "
             "problem. The operator code %d is not recognized", __func__,
@@ -981,6 +989,10 @@ reversepolish(struct arithmeticparams *p)
             { op=ARITHMETIC_OP_INTERPOLATE_MEDIANNGB; nop=0;  }
           else if (!strcmp(token->v, "collapse-sum"))
             { op=ARITHMETIC_OP_COLLAPSE_SUM;          nop=0; }
+          else if (!strcmp(token->v, "collapse-min"))
+            { op=ARITHMETIC_OP_COLLAPSE_MIN;          nop=0; }
+          else if (!strcmp(token->v, "collapse-max"))
+            { op=ARITHMETIC_OP_COLLAPSE_MAX;          nop=0; }
           else if (!strcmp(token->v, "collapse-mean"))
             { op=ARITHMETIC_OP_COLLAPSE_MEAN;         nop=0; }
           else if (!strcmp(token->v, "collapse-number"))
@@ -1081,6 +1093,8 @@ reversepolish(struct arithmeticparams *p)
                   break;
 
                 case ARITHMETIC_OP_COLLAPSE_SUM:
+                case ARITHMETIC_OP_COLLAPSE_MIN:
+                case ARITHMETIC_OP_COLLAPSE_MAX:
                 case ARITHMETIC_OP_COLLAPSE_MEAN:
                 case ARITHMETIC_OP_COLLAPSE_NUMBER:
                   arithmetic_collapse(p, token->v, op);
diff --git a/bin/arithmetic/arithmetic.h b/bin/arithmetic/arithmetic.h
index 67bff24..4aef527 100644
--- a/bin/arithmetic/arithmetic.h
+++ b/bin/arithmetic/arithmetic.h
@@ -42,6 +42,8 @@ enum arithmetic_prog_operators
   ARITHMETIC_OP_INVERT,
   ARITHMETIC_OP_INTERPOLATE_MEDIANNGB,
   ARITHMETIC_OP_COLLAPSE_SUM,
+  ARITHMETIC_OP_COLLAPSE_MIN,
+  ARITHMETIC_OP_COLLAPSE_MAX,
   ARITHMETIC_OP_COLLAPSE_MEAN,
   ARITHMETIC_OP_COLLAPSE_NUMBER,
 };
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 537ddb2..ad0bab2 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -11546,6 +11546,16 @@ values, all the elements in the returned dataset will 
have a single value
 (the length of the collapsed dimension). Therefore this is mostly relevant
 when there are blank values in the dataset.
 
address@hidden collapse-min
+Similar to @option{collapse-sum}, but the returned dataset will have the
+same numeric type as the input and will contain the minimum value for each
+pixel along the collapsed dimension.
+
address@hidden collapse-max
+Similar to @option{collapse-sum}, but the returned dataset will have the
+same numeric type as the input and will contain the maximum value for each
+pixel along the collapsed dimension.
+
 @item erode
 @cindex Erosion
 Erode the foreground pixels (with value @code{1}) of the input dataset
@@ -23231,6 +23241,15 @@ Arithmetic program's @option{collapse-number} operator 
(which uses this
 function) in @ref{Arithmetic operators}.
 @end deftypefun
 
address@hidden {gal_data_t *} gal_dimension_collapse_minmax (gal_data_t 
@code{*in}, size_t @code{c_dim}, int @code{max1_min0})
+Collapse the input dataset (@code{in}) along the given dimension
+(@code{c_dim}, in C definition: starting from zero, from the slowest
+dimension), by using the largest/smallest non-blank value along that
+dimension. If @code{max1_min0} is non-zero, then the collapsed dataset will
+have the maximum value along the given dimension and if it is zero, the
+minimum.
address@hidden deftypefun
+
 @deffn {Function-like macro} GAL_DIMENSION_NEIGHBOR_OP (@code{index}, 
@code{ndim}, @code{dsize}, @code{connectivity}, @code{dinc}, @code{operation})
 Parse the neighbors of the element located at @code{index} and do the
 requested operation on them. This is defined as a macro to allow easy
diff --git a/lib/dimension.c b/lib/dimension.c
index 9b1aea7..548b2c2 100644
--- a/lib/dimension.c
+++ b/lib/dimension.c
@@ -298,6 +298,8 @@ enum dimension_collapse_operation
  DIMENSION_COLLAPSE_INVALID,    /* ==0 by C standard. */
 
  DIMENSION_COLLAPSE_SUM,
+ DIMENSION_COLLAPSE_MAX,
+ DIMENSION_COLLAPSE_MIN,
  DIMENSION_COLLAPSE_MEAN,
  DIMENSION_COLLAPSE_NUMBER,
 };
@@ -371,22 +373,28 @@ dimension_collapse_sizes(gal_data_t *in, size_t c_dim, 
size_t *outndim,
 
 /* Depending on the operator, write the result into the output. */
 #define COLLAPSE_WRITE(OIND,IIND) {                                     \
-    /* We need the sum when number operator is requested. */            \
-    if(farr) farr[ OIND ] += (warr ? warr[w] : 1) * inarr[ IIND ];      \
+    /* Sum */                                                           \
+    if(farr)                                                            \
+      farr[ OIND ] += (warr ? warr[w] : 1) * inarr[ IIND ];             \
                                                                         \
-    /* We don't need the number in some situations. */                  \
+    /* Number */                                                        \
     if(iarr)                                                            \
       {                                                                 \
         if(num->type==GAL_TYPE_UINT8) iarr[ OIND ] = 1;                 \
         else                        ++iarr[ OIND ];                     \
       }                                                                 \
                                                                         \
-    /* If the sum of weights for is needed, add it. */                  \
+    /* Sum of weights. */                                               \
     if(wsumarr) wsumarr[ OIND ] += warr[w];                             \
+                                                                        \
+    /* Minimum or maximum. */                                           \
+    if(mmarr)                                                           \
+      mmarr[OIND] = ( max1_min0                                         \
+                      ? (mmarr[OIND]>inarr[IIND]?mmarr[OIND]:inarr[IIND]) \
+                      : (mmarr[OIND]<inarr[IIND]?mmarr[OIND]:inarr[IIND]) ); \
   }
 
 
-
 /* Deal properly with blanks. */
 #define COLLAPSE_CHECKBLANK(OIND,IIND) {                                \
     if(hasblank)                                                        \
@@ -406,8 +414,18 @@ dimension_collapse_sizes(gal_data_t *in, size_t c_dim, 
size_t *outndim,
 
 
 #define COLLAPSE_DIM(IT) {                                              \
-    IT B, *inarr=in->array;                                             \
+    IT m, B, *inarr=in->array, *mmarr=minmax?minmax->array:NULL;        \
     if(hasblank) gal_blank_write(&B, in->type);                         \
+                                                                        \
+    /* Initialize the array for minimum or maximum. */                  \
+    if(mmarr)                                                           \
+      {                                                                 \
+        if(max1_min0) gal_type_min(in->type, &m);                       \
+        else          gal_type_max(in->type, &m);                       \
+        for(i=0;i<minmax->size;++i) mmarr[i]=m;                         \
+      }                                                                 \
+                                                                        \
+    /* Collapse the dataset. */                                         \
     switch(in->ndim)                                                    \
       {                                                                 \
       /* 1D input dataset. */                                           \
@@ -458,6 +476,10 @@ dimension_collapse_sizes(gal_data_t *in, size_t c_dim, 
size_t *outndim,
               "supported, please contact us at %s to add this feature", \
               __func__, in->ndim, PACKAGE_BUGREPORT);                   \
       }                                                                 \
+                                                                        \
+    /* For minimum or maximum, elements with no input must be blank. */ \
+    if(mmarr && iarr)                                                   \
+      for(i=0;i<minmax->size;++i) if(iarr[i]==0) mmarr[i]=B;            \
   }
 
 
@@ -467,13 +489,14 @@ dimension_collapse_sizes(gal_data_t *in, size_t c_dim, 
size_t *outndim,
 gal_data_t *
 gal_dimension_collapse_sum(gal_data_t *in, size_t c_dim, gal_data_t *weight)
 {
+  int max1_min0=0;
   double *wsumarr=NULL;
   uint8_t *ii, *iarr=NULL;
   size_t a, b, i, j, k, w=-1, cnum=0;
   size_t outdsize[10], slice, outndim;
   int hasblank=gal_blank_present(in, 0);
   double *dd, *df, *warr=NULL, *farr=NULL;
-  gal_data_t *wht=NULL, *sum=NULL, *num=NULL;
+  gal_data_t *sum=NULL, *wht=NULL, *num=NULL, *minmax=NULL;
 
   /* Basic sanity checks. */
   wht=dimension_collapse_sanity_check(in, weight, c_dim, hasblank,
@@ -540,14 +563,15 @@ gal_data_t *
 gal_dimension_collapse_mean(gal_data_t *in, size_t c_dim,
                             gal_data_t *weight)
 {
+  int max1_min0=0;
   double wsum=NAN;
   double *wsumarr=NULL;
   int32_t *ii, *iarr=NULL;
   size_t a, b, i, j, k, w=-1, cnum=0;
   size_t outdsize[10], slice, outndim;
   int hasblank=gal_blank_present(in, 0);
-  gal_data_t *wht=NULL, *sum=NULL, *num=NULL;
   double *dd, *dw, *df, *warr=NULL, *farr=NULL;
+  gal_data_t *sum=NULL, *wht=NULL, *num=NULL, *minmax=NULL;
 
 
   /* Basic sanity checks. */
@@ -640,13 +664,14 @@ gal_dimension_collapse_mean(gal_data_t *in, size_t c_dim,
 gal_data_t *
 gal_dimension_collapse_number(gal_data_t *in, size_t c_dim)
 {
+  int max1_min0=0;
   double *wsumarr=NULL;
   double *warr=NULL, *farr=NULL;
   int32_t *ii, *iif, *iarr=NULL;
   size_t a, b, i, j, k, w, cnum=0;
   size_t outdsize[10], slice, outndim;
   int hasblank=gal_blank_present(in, 0);
-  gal_data_t *weight=NULL, *wht=NULL, *num=NULL;
+  gal_data_t *weight=NULL, *wht=NULL, *num=NULL, *minmax=NULL;
 
   /* Basic sanity checks. */
   wht=dimension_collapse_sanity_check(in, weight, c_dim, hasblank,
@@ -695,3 +720,65 @@ gal_dimension_collapse_number(gal_data_t *in, size_t c_dim)
   if(wht!=weight) gal_data_free(wht);
   return num;
 }
+
+
+
+
+
+gal_data_t *
+gal_dimension_collapse_minmax(gal_data_t *in, size_t c_dim, int max1_min0)
+{
+  int32_t *iarr=NULL;
+  double *wsumarr=NULL;
+  double *warr=NULL, *farr=NULL;
+  size_t a, b, i, j, k, w, cnum=0;
+  size_t outdsize[10], slice, outndim;
+  int hasblank=gal_blank_present(in, 0);
+  gal_data_t *weight=NULL, *wht=NULL, *num=NULL, *minmax=NULL;
+
+  /* Basic sanity checks. */
+  wht=dimension_collapse_sanity_check(in, weight, c_dim, hasblank,
+                                      &cnum, &warr);
+
+  /* Set the size of the collapsed output. */
+  dimension_collapse_sizes(in, c_dim, &outndim, outdsize);
+
+  /* Allocate the necessary datasets. If there are blank pixels, we'll need
+     to count how many elements whent into the calculation so we can set
+     them to blank. */
+  minmax=gal_data_alloc(NULL, in->type, outndim, outdsize, in->wcs,
+                        0, in->minmapsize, NULL, NULL, NULL);
+  if(hasblank)
+    {
+      num=gal_data_alloc(NULL, GAL_TYPE_INT32, outndim, outdsize, in->wcs,
+                         1, in->minmapsize, NULL, NULL, NULL);
+      iarr=num->array;
+    }
+
+  /* Parse the input dataset (if necessary). */
+  switch(in->type)
+    {
+    case GAL_TYPE_UINT8:     COLLAPSE_DIM( uint8_t  );   break;
+    case GAL_TYPE_INT8:      COLLAPSE_DIM( int8_t   );   break;
+    case GAL_TYPE_UINT16:    COLLAPSE_DIM( uint16_t );   break;
+    case GAL_TYPE_INT16:     COLLAPSE_DIM( int16_t  );   break;
+    case GAL_TYPE_UINT32:    COLLAPSE_DIM( uint32_t );   break;
+    case GAL_TYPE_INT32:     COLLAPSE_DIM( int32_t  );   break;
+    case GAL_TYPE_UINT64:    COLLAPSE_DIM( uint64_t );   break;
+    case GAL_TYPE_INT64:     COLLAPSE_DIM( int64_t  );   break;
+    case GAL_TYPE_FLOAT32:   COLLAPSE_DIM( float    );   break;
+    case GAL_TYPE_FLOAT64:   COLLAPSE_DIM( double   );   break;
+    default:
+      error(EXIT_FAILURE, 0, "%s: type value (%d) not recognized",
+            __func__, in->type);
+    }
+
+  /* Remove the respective dimension in the WCS structure also (if any
+     exists). Note that `sum->ndim' has already been changed. So we'll use
+     `in->wcs'. */
+  gal_wcs_remove_dimension(minmax->wcs, in->ndim-c_dim);
+
+  /* Return. */
+  if(wht!=weight) gal_data_free(wht);
+  return minmax;
+}
diff --git a/lib/gnuastro/dimension.h b/lib/gnuastro/dimension.h
index a0623a2..46e68fa 100644
--- a/lib/gnuastro/dimension.h
+++ b/lib/gnuastro/dimension.h
@@ -109,7 +109,8 @@ gal_dimension_collapse_mean(gal_data_t *in, size_t c_dim,
 gal_data_t *
 gal_dimension_collapse_number(gal_data_t *in, size_t c_dim);
 
-
+gal_data_t *
+gal_dimension_collapse_minmax(gal_data_t *in, size_t c_dim, int max1_min0);