[gnuastro-commits] master 465eac4e 2/3: Arithmetic: new operator called add-dimension-fast

[Mohammad Akhlaghi]

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

    Arithmetic: new operator called add-dimension-fast
    
    Until now, there was only one operator to add multiple datasets into a
    higher-dimensional dataset: 'add-dimension'. But it would only work on 2D
    images, and add datasets along the slowest dimension.
    
    With this commit, a new 'add-dimension-fast' operator has been added that
    currently only takes 1D datasets, but can add new datasets along the
    fastest dimension. To avoid confusion with the existing function, the
    existing function has been renamed to 'add-dimension-slow'.
    
    In the process of testing this new feature, I also recognized that in
    'gal_type_string_to_number', we were mistakenly using 'abs' (for integers)
    to check the range of 32-bit floating point values. It has been corrected
    to 'fabs' (for double type).
---
 NEWS                        |  11 ++++
 bin/arithmetic/arithmetic.c | 119 ++++++++++++++++++++++++++++++++------------
 bin/arithmetic/arithmetic.h |   3 +-
 doc/gnuastro.texi           |  14 ++++--
 lib/type.c                  |   2 +-
 5 files changed, 111 insertions(+), 38 deletions(-)

diff --git a/NEWS b/NEWS
index 0d469609..f14cd1bc 100644
--- a/NEWS
+++ b/NEWS
@@ -31,6 +31,14 @@ See the end of the file for license conditions.
      improve the running time of the programs. This task was implemented
      based on feedback from Andrés Del Pino Molina.
 
+  Arithmetic:
+   - 'add-dimension-fast': Add a new dataset along the "fastest" dimension
+     of the first dataset (in a FITS image, the "fast" dimension is the
+     horizontal one). For example if you have N one-dimensional datasets
+     with M elements, you can use this operator to construct a
+     2-dimensional FITS image that has N pixels along the horizontal and M
+     pixels along the vertical.
+
   Crop:
    --widthinpix: when in WCS-mode, the value given to '--width' will be
      interpreted in units of pixels. This will allow having crops of a
@@ -92,6 +100,9 @@ See the end of the file for license conditions.
 
 ** Changed features
 
+  Arithmetic:
+   - 'add-dimension-slow' is the new name for the old 'add-dimension'.
+
   Segment:
    - '--noobjects' is the new name for the old '--onlyclumps' option. This
      is because the name Only-Clumps can lead to a wrong expectation that
diff --git a/bin/arithmetic/arithmetic.c b/bin/arithmetic/arithmetic.c
index ca85023a..6d1bc1c9 100644
--- a/bin/arithmetic/arithmetic.c
+++ b/bin/arithmetic/arithmetic.c
@@ -1056,18 +1056,20 @@ arithmetic_unique(struct arithmeticparams *p, char 
*token, int operator)
 
 
 void
-arithmetic_add_dimension(struct arithmeticparams *p, char *token, int operator)
+arithmetic_add_dimension(struct arithmeticparams *p, char *token,
+                         int operator)
 {
-  gal_data_t *out=NULL;
+  size_t one=1;
+  gal_data_t *ref=NULL, *out=NULL;
+  size_t i, j, num, dsize[3], nbytes=0;
   gal_data_t *tmp = operands_pop(p, token);
-  size_t i, num, dsize[3], ndim=3, nbytes=0;
 
   /* Make sure the first operand is a number. */
   if(tmp->size!=1)
     error(EXIT_FAILURE, 0, "first popped operand to '%s' must be a "
           "number (specifying how many datasets to use)", token);
 
-  /* Put the value into 'num'. */
+  /* Put the first-popped value into 'num'. */
   tmp=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_SIZE_T);
   num=*(size_t *)(tmp->array);
   gal_data_free(tmp);
@@ -1081,43 +1083,88 @@ arithmetic_add_dimension(struct arithmeticparams *p, 
char *token, int operator)
       /* Things that differ from the first dataset and the rest. */
       if(out) /* Not the first. */
         {
-          /* Basic sanity checks. */
+          /* All entries have to have the same data type. */
           if(tmp->type!=out->type)
-            error(EXIT_FAILURE, 0, "the operands to '%s' have to have the "
-                  "same data type (the inputs contain atleast two types: "
-                  "'%s' and '%s')", token, gal_type_name(tmp->type, 1),
-                  gal_type_name(out->type, 1));
-          if( tmp->ndim!=out->ndim-1
-              || tmp->dsize[0]!=out->dsize[1]
-              || tmp->dsize[1]!=out->dsize[2] )
-            error(EXIT_FAILURE, 0, "the operands to '%s' have to have the "
-                  "same size", token);
+            error(EXIT_FAILURE, 0, "the operands to '%s' must have "
+                  "the same data type (the first popped operand has a "
+                  "'%s' type, but the popped operand number %zu has "
+                  "type '%s')", token, gal_type_name(out->type, 1), i+1,
+                  gal_type_name(tmp->type, 1));
+
+          /* All entries should also have the same dimension. */
+          if( gal_dimension_is_different(ref, tmp) )
+            error(EXIT_FAILURE, 0, "the operands to '%s' must have the "
+                  "same size in all dimensions", token);
         }
       else  /* First popped operand. */
         {
           /* First popped operand, do necessary basic checks here. */
-          if(tmp->ndim!=2)
-            error(EXIT_FAILURE, 0, "currently only 2-dimensional datasets "
-                  "are acceptable for '%s', please get in touch with us at "
-                  "%s so we add functionality for different dimensions",
-                  token, PACKAGE_BUGREPORT);
+          switch(operator)
+            {
+            case ARITHMETIC_OP_ADD_DIMENSION_SLOW:
+              dsize[0]=num;
+              dsize[1]=tmp->dsize[0];
+              dsize[2]=tmp->dsize[1];
+              nbytes=gal_type_sizeof(tmp->type)*tmp->size;
+              if(tmp->ndim!=2)
+                error(EXIT_FAILURE, 0, "currently only 2-dimensional "
+                      "datasets are acceptable for '%s', please get in "
+                      "touch with us at %s so we add functionality for "
+                      "different dimensions", token, PACKAGE_BUGREPORT);
+              break;
+            case ARITHMETIC_OP_ADD_DIMENSION_FAST:
+              dsize[1]=num;
+              dsize[0]=tmp->dsize[0];
+              nbytes=gal_type_sizeof(tmp->type);
+              if(tmp->ndim!=1)
+                error(EXIT_FAILURE, 0, "currently only 1-dimensional "
+                      "datasets are acceptable for '%s', please get in "
+                      "touch with us at %s so we add functionality for "
+                      "different dimensions", token, PACKAGE_BUGREPORT);
+              break;
+            default:
+              error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at "
+                    "%s to address the problem. The operator code '%d' "
+                    "isn't recognized", __func__, PACKAGE_BUGREPORT,
+                    operator);
+            }
+
+          /* Allocate the size-reference dataset to simplify checking of
+             the sizes for each new dataset. This is done because this
+             dataset that is actually being popped will be freed right
+             after its values have been written into the output. We don't
+             care about any data from this dataset, and to avoid wasting
+             memory, we will simply allocate a one-element dataset. But
+             we'll many set its sizes to the size of the popped dataset. */
+          ref=gal_data_alloc(NULL, tmp->type, 1, &one, NULL, 0, -1, 1,
+                             NULL, NULL, NULL);
+          free(ref->dsize);
+          ref->ndim=tmp->ndim;
+          ref->dsize=gal_pointer_allocate(GAL_TYPE_SIZE_T, ref->ndim, 0,
+                                          __func__, "ref->dsize");
+          for(j=0;j<ref->ndim;++j) ref->dsize[j]=tmp->dsize[j];
 
           /* Allocate the output dataset. */
-          dsize[0]=num;
-          dsize[1]=tmp->dsize[0];
-          dsize[2]=tmp->dsize[1];
-          out = gal_data_alloc(NULL, tmp->type, ndim, dsize, NULL, 0,
-                               p->cp.minmapsize, p->cp.quietmmap, NULL,
+          out = gal_data_alloc(NULL, tmp->type, tmp->ndim+1, dsize, NULL,
+                               0, p->cp.minmapsize, p->cp.quietmmap, NULL,
                                NULL, NULL);
-
-          /* Get the number of bytes in each dataset. */
-          nbytes=gal_type_sizeof(tmp->type)*tmp->size;
         }
 
       /* Copy the dataset into the higher-dimensional output. */
-      memcpy(gal_pointer_increment(out->array, (num-1-i)*tmp->size,
-                                   tmp->type),
-             tmp->array, nbytes);
+      switch(operator)
+        {
+        case ARITHMETIC_OP_ADD_DIMENSION_SLOW:
+          memcpy(gal_pointer_increment(out->array, (num-1-i)*tmp->size,
+                                       tmp->type),
+                 tmp->array, nbytes);
+          break;
+        case ARITHMETIC_OP_ADD_DIMENSION_FAST:
+          for(j=0;j<tmp->size;++j)
+            memcpy(gal_pointer_increment(out->array, j*num+num-i-1, out->type),
+                   gal_pointer_increment(tmp->array, j, out->type), nbytes);
+
+          break;
+        }
 
       /* Clean up. */
       gal_data_free(tmp);
@@ -1125,6 +1172,9 @@ arithmetic_add_dimension(struct arithmeticparams *p, char 
*token, int operator)
 
   /* Put the higher-dimensional output on the operands stack. */
   operands_add(p, NULL, out);
+
+  /* Clean up. */
+  gal_data_free(ref);
 }
 
 
@@ -1199,8 +1249,10 @@ arithmetic_set_operator(char *string, size_t 
*num_operands)
         { op=ARITHMETIC_OP_COLLAPSE_NUMBER;       *num_operands=0; }
       else if (!strcmp(string, "unique"))
         { op=ARITHMETIC_OP_UNIQUE;                *num_operands=0; }
-      else if (!strcmp(string, "add-dimension"))
-        { op=ARITHMETIC_OP_ADD_DIMENSION;         *num_operands=0; }
+      else if (!strcmp(string, "add-dimension-slow"))
+        { op=ARITHMETIC_OP_ADD_DIMENSION_SLOW;    *num_operands=0; }
+      else if (!strcmp(string, "add-dimension-fast"))
+        { op=ARITHMETIC_OP_ADD_DIMENSION_FAST;    *num_operands=0; }
       else
         error(EXIT_FAILURE, 0, "the argument '%s' could not be "
               "interpretted as a file name, named dataset, number, "
@@ -1334,7 +1386,8 @@ arithmetic_operator_run(struct arithmeticparams *p, int 
operator,
           arithmetic_unique(p, operator_string, operator);
           break;
 
-        case ARITHMETIC_OP_ADD_DIMENSION:
+        case ARITHMETIC_OP_ADD_DIMENSION_SLOW:
+        case ARITHMETIC_OP_ADD_DIMENSION_FAST:
           arithmetic_add_dimension(p, operator_string, operator);
           break;
 
diff --git a/bin/arithmetic/arithmetic.h b/bin/arithmetic/arithmetic.h
index 4ed95a8d..62828c56 100644
--- a/bin/arithmetic/arithmetic.h
+++ b/bin/arithmetic/arithmetic.h
@@ -51,7 +51,8 @@ enum arithmetic_prog_operators
   ARITHMETIC_OP_COLLAPSE_MEAN,
   ARITHMETIC_OP_COLLAPSE_NUMBER,
   ARITHMETIC_OP_UNIQUE,
-  ARITHMETIC_OP_ADD_DIMENSION,
+  ARITHMETIC_OP_ADD_DIMENSION_SLOW,
+  ARITHMETIC_OP_ADD_DIMENSION_FAST,
 };
 
 
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index a8ecb5da..6f346ca9 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -14451,7 +14451,7 @@ Similar to @option{collapse-sum}, but the returned 
dataset will have the same nu
 @item 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 add-dimension
+@item add-dimension-slow
 Build a higher-dimensional dataset from all the input datasets stacked after 
one another (along the slowest dimension).
 The first popped operand has to be a single number.
 It is used by the operator to know how many operands it should pop from the 
stack (and the size of the output in the new dimension).
@@ -14464,11 +14464,19 @@ If no file is specified for the WCS, the first 
dataset's WCS will be used, you c
 If your datasets don't have the same type, you can use the type transformation 
operators of Arithmetic that are discussed below.
 Just beware of overflow if you are transforming to a smaller type, see 
@ref{Numeric data types}.
 
-For example if you want to put the three @file{img1.fits}, @file{img2.fits} 
and @file{img3.fits} images (each a 2D dataset) into one 3D datacube, you can 
use this command:
+For example, let's assume you have 3 two-dimensional images @file{a.fits}, 
@file{b.fits} and @file{c.fits} (each with @mymath{200\times100} pixels).
+You can construct a 3D data cube with @mymath{200\times100\times3} voxels 
(volume-pixels) using the command below:
 
 @example
-$ astarithmetic img1.fits img2.fits img3.fits 3 add-dimension
+$ astarithmetic a.fits b.fits c.fits 3 add-dimension-slow
 @end example
+
+@item add-dimension-fast
+Similar to @code{add-dimension-slow} but along the fastest dimension.
+This operator currently only works for 1D input operands, please contact us if 
you want inputs to have different dimensions.
+
+For example, let's assume you have 3 one-dimensional datasets, each with 100 
elements.
+With this operator, you can construct a @mymath{3\times100} pixel FITS image 
that has 3 pixels along the horizontal and 5 pixels along the vertical.
 @end table
 
 @node Conditional operators, Mathematical morphology operators, Dimensionality 
changing operators, Arithmetic operators
diff --git a/lib/type.c b/lib/type.c
index a6dd8ba0..89e8e501 100644
--- a/lib/type.c
+++ b/lib/type.c
@@ -639,7 +639,7 @@ gal_type_string_to_number(char *string, uint8_t *type)
 
       /* Calculate the number of decimal digits and decide if it the number
          should be a float or a double. */
-      if( lnz-fnz > FLT_DIG || abs(d)>FLT_MAX || abs(d)<FLT_MIN )
+      if( lnz-fnz > FLT_DIG || fabs(d)>FLT_MAX || fabs(d)<FLT_MIN )
         {      ptr=&d; *type=GAL_TYPE_FLOAT64; }
       else
         { f=d; ptr=&f; *type=GAL_TYPE_FLOAT32; }