[gnuastro-commits] master ddcb77a 2/2: Arithmetic library: more complete sanity checks in arithmetic_size

[Mohammad Akhlaghi]

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

    Arithmetic library: more complete sanity checks in arithmetic_size
    
    Until now, the `arithmetic_size' function had the following issues:
      1) The sanity checks weren't complete.
      2) the `arg' dataset wouldn't be freed/not-freed properly.
      3) it would report the length along the C dimension, while astronomers
         expect the FITS order (which is similar to Fortran's dimension
         ordering (the slowest dimension is counted as first).
    
    With this commit, some extra sanity checks were added to the
    `arithmetic_size' function of `lib/arithmetic.c', and the error messages of
    the existing ones were made more clear for the users. Also, the value of
    the `arg' dataset is now immediately written into a `size_t' variable, and
    the dataset is freed (or not!) at that point. Finally the output value is
    the length along the FITS order.
    
    Also in the documentation, the description of the new `size' operator was
    moved after all the operators that deal with a dataset's data. This is
    because this operator deals with the dataset's metadata. Later, other
    operators that deal with metadata should be added under it. We should also
    generally start classifying the operators to be more readable/findable.
---
 NEWS              |  1 +
 doc/gnuastro.texi | 19 +++++++-------
 lib/arithmetic.c  | 75 ++++++++++++++++++++++++++++++++-----------------------
 3 files changed, 55 insertions(+), 40 deletions(-)

diff --git a/NEWS b/NEWS
index 0704273..b6e84ac 100644
--- a/NEWS
+++ b/NEWS
@@ -9,6 +9,7 @@ See the end of the file for license conditions.
 
   Arithmetic:
    - New `quantile' operator for coadding datasets.
+   - New `size' operator to report length of dataset in requested dimension.
 
   CosmicCalculator:
    --listlines: list the pre-defined spectral line wavelengths and names
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 53b2a9e..e6414dd 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -10154,15 +10154,6 @@ 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 size
-Size of the dataset along a given dimension. 
-The output will be a single unsigned integer(dimensions cannot be negative).
-For example, the following command will produce the dimension of the 
hdu1(default) in the first operand(a.fits) along the second axis.
-
-@example
-astarithmetic a.fits 2 size
-@end example
-
 @item filter-mean
 Apply mean filtering (or @url{https://en.wikipedia.org/wiki/Moving_average, 
moving average}) on the input dataset.
 During mean filtering, each pixel (data element) is replaced by the mean value 
of all its surrounding pixels (excluding blank values).
@@ -10497,6 +10488,16 @@ The internal conversion of C will be used.
 Convert the type of the popped operand to 64-bit (double precision) floating 
point (see @ref{Numeric data types}).
 The internal conversion of C will be used.
 
+@item size
+Size of the dataset along a given FITS/Fortran dimension (counting from 1).
+The desired dimension should be the first popped operand and the dataset must 
be the second popped operand.
+The output will be a single unsigned integer (dimensions cannot be negative).
+For example, the following command will produce the size of the first 
extension/HDU (the default HDU) of @file{a.fits} along the second FITS axis.
+
+@example
+astarithmetic a.fits 2 size
+@end example
+
 @item set-AAA
 Set the characters after the dash (@code{AAA} in the case shown here) as a 
name for the first popped operand on the stack.
 The named dataset will be freed from memory as soon as it is no longer needed, 
or if the name is reset to refer to another dataset later in the command.
diff --git a/lib/arithmetic.c b/lib/arithmetic.c
index b759428..bab5707 100644
--- a/lib/arithmetic.c
+++ b/lib/arithmetic.c
@@ -558,47 +558,60 @@ arithmetic_from_statistics(int operator, int flags, 
gal_data_t *input)
 static gal_data_t *
 arithmetic_size(int operator, int flags, gal_data_t *in, gal_data_t *arg)
 {
-  gal_data_t *out=NULL;
-  size_t *out_arr;
-  size_t arg_val, temp_dsize=1;
-  size_t *dsize=&temp_dsize;
-
-  /* Check the type */
-  switch(arg->type)
+  size_t one=1, arg_val;
+  gal_data_t *usearg=NULL, *out=NULL;
+
+  /* Sanity checks on argument (dimension number): it should be an integer,
+     and have a size of 1. */
+  if(arg->type==GAL_TYPE_FLOAT32 || arg->type==GAL_TYPE_FLOAT64)
+    error(EXIT_FAILURE, 0, "%s: size operator's dimention argument"
+          "must have an interger type", __func__);
+  if(arg->size!=1)
+    error(EXIT_FAILURE, 0, "%s: size operator's dimention argument"
+          "must be a single number, but it has %zu elements", __func__,
+          arg->size);
+
+
+  /* Convert `arg' to `size_t' and read it. Note that we can only free the
+     `arg' array (while changing its type), when the freeing flag has been
+     set. */
+  if(flags & GAL_ARITHMETIC_FREE)
     {
-      case GAL_TYPE_FLOAT32:
-      case GAL_TYPE_FLOAT64:
-        error(EXIT_FAILURE, 0, "%s: size "
-              "operator can only work on integer types", __func__);
+      arg=gal_data_copy_to_new_type_free(arg, GAL_TYPE_SIZE_T);
+      arg_val=*(size_t *)(arg->array);
+      gal_data_free(arg);
+    }
+  else
+    {
+      usearg=gal_data_copy_to_new_type(arg, GAL_TYPE_SIZE_T);
+      arg_val=*(size_t *)(usearg->array);
+      gal_data_free(usearg);
     }
 
-  /* Convert `arg` to GAL_TYPE_SIZE_T */
-  arg=gal_data_copy_to_new_type_free(arg, GAL_TYPE_SIZE_T);
 
-  arg_val=*(size_t *)(arg->array);
+  /* Sanity checks on the value of the given argument.*/
+  if(arg_val>in->ndim)
+    error(EXIT_FAILURE, 0, "%s: size operator's dimension argument "
+          "(given %zu) cannot be larger than the dimentions of the "
+          "given input (%zu)", __func__, arg_val, in->ndim);
+  if(arg_val==0)
+    error(EXIT_FAILURE, 0, "%s: size operator's dimension argument "
+          "(given %zu) cannot be zero: dimensions are counted from 1",
+          __func__, arg_val);
 
-  /* Allocate size for the output data to make changes.
-     Set the output array value. */
-  dsize[0]=temp_dsize;
-  out=gal_data_alloc(NULL, GAL_TYPE_UINT64, 1, dsize,
-                    NULL, 0, in->minmapsize, 0,
-                    NULL, NULL, NULL);
 
-  out_arr=out->array;
-  out_arr[0]=in->dsize[arg_val-1];
+  /* Allocate the output array and write the desired dimension. Note that
+     `dsize' is in the C order, while the output must be in FITS/Fortran
+     order. Also that C order starts from 0, while the FITS order starts
+     from 1. */
+  out=gal_data_alloc(NULL, GAL_TYPE_SIZE_T, 1, &one, NULL, 0,
+                     in->minmapsize, 0, NULL, NULL, NULL);
+  *(size_t *)(out->array)=in->dsize[in->ndim-arg_val];
 
-  /* The first argument must be a single integer. */
-  if(arg->size != 1)
-    {
-      /* If an image/array is passed, report error. */
-      error(EXIT_FAILURE, 0, "%s: first argument must be a single integer. ",
-            __func__);
-    }
 
   /* Clean up and return */
-  if( flags & GAL_ARITHMETIC_FREE)
+  if(flags & GAL_ARITHMETIC_FREE)
     gal_data_free(in);
-
   return out;
 }