[gnuastro-commits] master 09359bc: Arithmetic: add-dimension operator to build higher-dimensional outputs

[Mohammad Akhlaghi]

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

    Arithmetic: add-dimension operator to build higher-dimensional outputs
    
    With this operator, it is possible to put individual datasets into a
    higher-dimensional dataset.
---
 NEWS                        |  3 ++
 bin/arithmetic/arithmetic.c | 84 ++++++++++++++++++++++++++++++++++++++++++++-
 bin/arithmetic/arithmetic.h |  1 +
 doc/gnuastro.texi           | 19 ++++++++++
 4 files changed, 106 insertions(+), 1 deletion(-)

diff --git a/NEWS b/NEWS
index 5cca4f3..7cb6478 100644
--- a/NEWS
+++ b/NEWS
@@ -8,6 +8,9 @@ See the end of the file for license conditions.
 ** New features
 
   Arithmetic:
+   - The new `add-dimension' operator will stack the popped operands in a
+     higher-dimensional dataset. For example to build a 3D cube from
+     individual 2D images/slices.
    --onedonstdout: when the output is one-dimensional, print the values on
      the standard output, not into a file.
 
diff --git a/bin/arithmetic/arithmetic.c b/bin/arithmetic/arithmetic.c
index 6935036..82bbba3 100644
--- a/bin/arithmetic/arithmetic.c
+++ b/bin/arithmetic/arithmetic.c
@@ -928,7 +928,7 @@ arithmetic_unique(struct arithmeticparams *p, char *token, 
int operator)
   /* Remove all blank elements. */
   gal_blank_remove(input);
 
-  /* Clean up and add the collapsed dataset to the top of the operands. */
+  /* Add the collapsed dataset to the top of the operands. */
   operands_add(p, NULL, input);
 }
 
@@ -936,6 +936,82 @@ arithmetic_unique(struct arithmeticparams *p, char *token, 
int operator)
 
 
 
+void
+arithmetic_add_dimension(struct arithmeticparams *p, char *token, int operator)
+{
+  gal_data_t *out=NULL;
+  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'. */
+  tmp=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_SIZE_T);
+  num=*(size_t *)(tmp->array);
+  gal_data_free(tmp);
+
+  /* Pop all the datasets and put them in a list. */
+  for(i=0;i<num;++i)
+    {
+      /* Pop the operand. */
+      tmp=operands_pop(p, token);
+
+      /* Things that differ from the first dataset and the rest. */
+      if(out) /* Not the first. */
+        {
+          /* Basic sanity checks. */
+          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);
+        }
+      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);
+
+          /* 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,
+                               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);
+
+      /* Clean up. */
+      gal_data_free(tmp);
+    }
+
+  /* Put the higher-dimensional output on the operands stack. */
+  operands_add(p, NULL, out);
+}
+
+
+
+
+
 
 
 
@@ -996,6 +1072,8 @@ 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
         error(EXIT_FAILURE, 0, "the argument \"%s\" could not be "
               "interpretted as a file name, named dataset, number, "
@@ -1119,6 +1197,10 @@ arithmetic_operator_run(struct arithmeticparams *p, int 
operator,
           arithmetic_unique(p, operator_string, operator);
           break;
 
+        case ARITHMETIC_OP_ADD_DIMENSION:
+          arithmetic_add_dimension(p, operator_string, operator);
+          break;
+
         default:
           error(EXIT_FAILURE, 0, "%s: a bug! please contact us at "
                 "%s to fix the problem. The code %d is not "
diff --git a/bin/arithmetic/arithmetic.h b/bin/arithmetic/arithmetic.h
index d267a58..4e90d99 100644
--- a/bin/arithmetic/arithmetic.h
+++ b/bin/arithmetic/arithmetic.h
@@ -47,6 +47,7 @@ enum arithmetic_prog_operators
   ARITHMETIC_OP_COLLAPSE_MEAN,
   ARITHMETIC_OP_COLLAPSE_NUMBER,
   ARITHMETIC_OP_UNIQUE,
+  ARITHMETIC_OP_ADD_DIMENSION,
 };
 
 
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 693ed79..78c8b22 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -9905,6 +9905,25 @@ 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
+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).
+The rest of the operands must have the same size and numerical data type.
+This operator currently only works for 2D input operands, please contact us if 
you want inputs to have different dimensions.
+
+The output's WCS (which should have a different dimensionality compared to the 
inputs) can be read from another file with the @option{--wcsfile} option.
+If no file is specified for the WCS, the first dataset's WCS will be used, you 
can later add/change the necessary WCS keywords with the FITS keyword 
modification features of the Fits program (see @ref{Fits}).
+
+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:
+
+@example
+$ astarithmetic img1.fits img2.fits img3.fits 3 add-dimension
+@end example
+
 @item unique
 Remove all duplicate (and blank) elements from the first popped operand.
 The unique elements of the dataset will be stored in a single-dimensional 
dataset.