[gnuastro-commits] master 44e27da2: Library (arithmetic): operators properly account for empty datasets

[Mohammad Akhlaghi]

branch: master
commit 44e27da23669fc77eb58f4b9320347f6f0a59f5b
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>

    Library (arithmetic): operators properly account for empty datasets
    
    Until now, all the arithmetic functions in the library (and the
    column-specific operators of Table's column arithmetic!), implicitly
    assumed that an actual array exists in the input dataset. However,
    especially when dealing with Column-arithmetic in Table, this is not always
    the case: it can happen that the input table is from a previous automatic
    command where no rows were selected!
    
    In such cases, Table would abort with a "Segmentation fault (core dumped)",
    which is not good and not informative!
    
    With this commit a check has been added at the start of all the arithmetic
    library and column-specific functions that will properly deal with empty
    datasets.
    
    This bug was reported by Zohreh Ghaffari.
    
    This fixes bug #62204.
---
 NEWS                   |  2 ++
 bin/table/arithmetic.c | 61 ++++++++++++++++++++---------------
 lib/arithmetic.c       | 87 ++++++++++++++++++++++++++++++++++++++++++++++++--
 lib/blank.c            |  5 +++
 lib/data.c             | 17 ++++++++--
 lib/wcs.c              | 17 ++++++++--
 6 files changed, 155 insertions(+), 34 deletions(-)

diff --git a/NEWS b/NEWS
index 137159ef..c0ad9d36 100644
--- a/NEWS
+++ b/NEWS
@@ -263,6 +263,8 @@ See the end of the file for license conditions.
               (that guarantees monotonicity) not used in Gnuastro's
               library, even if it is present and found at configure time;
               found and fixed by Pedram Ashofteh Ardakani.
+  bug #62204: Column arithmetic crashes when there are no rows; reported by
+              Zohreh Ghaffari.
 
 
 
diff --git a/bin/table/arithmetic.c b/bin/table/arithmetic.c
index e2af7753..4e88c841 100644
--- a/bin/table/arithmetic.c
+++ b/bin/table/arithmetic.c
@@ -622,17 +622,20 @@ arithmetic_distance(struct tableparams *p, gal_data_t 
**stack, int operator)
                      p->cp.minmapsize, p->cp.quietmmap, colname, NULL,
                      colcomment);
 
-  /* Measure the distances.  */
-  o=out->array;
-  a1=a->array; a2=a->next->array;
-  b1=b->array; b2=b->next->array;
-  if(a->size==1 || b->size==1) /* One of them is a single point. */
-    for(i=0;i<a->size;++i)
-      for(j=0;j<b->size;++j)
-        o[a->size>b->size?i:j] = distance_func(a1[i], a2[i], b1[j], b2[j]);
-  else                         /* Both have the same length. */
-    for(i=0;i<a->size;++i)     /* (all were originally from the same table) */
-      o[i] = distance_func(a1[i], a2[i], b1[i], b2[i]);
+  /* Measure the distances (if the dataset isn't empty: it can be!). */
+  if(out->size>0 && out->array)
+    {
+      o=out->array;
+      a1=a->array; a2=a->next->array;
+      b1=b->array; b2=b->next->array;
+      if(a->size==1 || b->size==1) /* One of them is a single point. */
+        for(i=0;i<a->size;++i)
+          for(j=0;j<b->size;++j)
+            o[a->size>b->size?i:j] = distance_func(a1[i], a2[i], b1[j], b2[j]);
+      else                     /* Both have the same length. */
+        for(i=0;i<a->size;++i) /* (all were originally from the same table) */
+          o[i] = distance_func(a1[i], a2[i], b1[i], b2[i]);
+    }
 
   /* Clean up and put the output dataset onto the stack. */
   gal_list_data_free(a);
@@ -690,21 +693,26 @@ arithmetic_datetosec(struct tableparams *p, gal_data_t 
**stack,
                      p->cp.minmapsize, p->cp.quietmmap, name, unit,
                      comment);
 
-  /* Convert each input string into number of seconds. */
-  iarr=out->array;
-  for(i=0; i<in->size; ++i)
+  /* Convert each input string into number of seconds. Note that it is
+     possible to have an empty dataset, in that case, we shouldn't do
+     anything.*/
+  if(out->size>0 && out->array)
     {
-      /* Read the number of seconds and sub-seconds and write into the
-         output. */
-      v=gal_fits_key_date_to_seconds(strarr[i], &subsecstr,
-                                     &subsec);
-      iarr[i] = ( v==GAL_BLANK_SIZE_T
-                  ? GAL_BLANK_INT64
-                  : ( operator == ARITHMETIC_TABLE_OP_DATETOSEC
-                      ? v
-                      : (isnan(subsec)
-                         ? v*1000
-                         : v*1000 + (int64_t)(subsec*1000) ) ) );
+      iarr=out->array;
+      for(i=0; i<in->size; ++i)
+        {
+          /* Read the number of seconds and sub-seconds and write into the
+             output. */
+          v=gal_fits_key_date_to_seconds(strarr[i], &subsecstr,
+                                         &subsec);
+          iarr[i] = ( v==GAL_BLANK_SIZE_T
+                      ? GAL_BLANK_INT64
+                      : ( operator == ARITHMETIC_TABLE_OP_DATETOSEC
+                          ? v
+                          : (isnan(subsec)
+                             ? v*1000
+                             : v*1000 + (int64_t)(subsec*1000) ) ) );
+        }
     }
 
   /* Clean up and put the resulting calculation back on the stack. */
@@ -819,7 +827,8 @@ arithmetic_operator_run(struct tableparams *p,
          arguments it uses depend on the operator. In other words, when the
          operator doesn't need three operands, the extra arguments will be
          ignored. */
-      gal_list_data_add(stack, gal_arithmetic(token->operator, 
p->cp.numthreads,
+      gal_list_data_add(stack, gal_arithmetic(token->operator,
+                                              p->cp.numthreads,
                                               flags, d1, d2, d3) );
 
       /* Reset the meta-data for the element that was just put on the
diff --git a/lib/arithmetic.c b/lib/arithmetic.c
index 4eac9746..6f09ba83 100644
--- a/lib/arithmetic.c
+++ b/lib/arithmetic.c
@@ -182,6 +182,11 @@ arithmetic_not(gal_data_t *data, int flags)
   uint8_t *o;
   gal_data_t *out;
 
+  /* The dataset may be empty. In this case the output should also be empty
+     (we can have tables and images with 0 rows or pixels!). */
+  if(data->size==0 || data->array==NULL) return data;
+
+
   /* Allocate the output array. */
   out=gal_data_alloc(NULL, GAL_TYPE_UINT8, data->ndim, data->dsize,
                      data->wcs, 0, data->minmapsize, data->quietmmap,
@@ -238,6 +243,10 @@ arithmetic_bitwise_not(int flags, gal_data_t *in)
   uint64_t   *iu64 = in->array,  *iu64f = iu64 + in->size,   *ou64;
   int64_t    *ii64 = in->array,  *ii64f = ii64 + in->size,   *oi64;
 
+  /* The dataset may be empty. In this case, the output should also be
+     empty (we can have tables and images with 0 rows or pixels!). */
+  if(in->size==0 || in->array==NULL) return in;
+
   /* Check the type */
   switch(in->type)
     {
@@ -313,6 +322,10 @@ arithmetic_abs(int flags, gal_data_t *in)
 {
   gal_data_t *out;
 
+  /* The dataset may be empty. In this case, the output should also be
+     empty (we can have tables and images with 0 rows or pixels!). */
+  if(in->size==0 || in->array==NULL) return in;
+
   /* Set the output array. */
   if(flags & GAL_ARITHMETIC_FLAG_INPLACE)
     out=in;
@@ -502,6 +515,10 @@ arithmetic_function_unary(int operator, int flags, 
gal_data_t *in)
   gal_data_t *o;
   double pi=3.14159265358979323846264338327;
 
+  /* The dataset may be empty. In this case, the output should also be empty
+     (we can have tables and images with 0 rows or pixels!). */
+  if(in->size==0 || in->array==NULL) return in;
+
   /* See if the operation should be done in place. The output of these
      operators is defined in the floating point space. So even if the input
      is integer type and user requested inplace opereation, if its not a
@@ -626,8 +643,14 @@ static gal_data_t *
 arithmetic_from_statistics(int operator, int flags, gal_data_t *input)
 {
   gal_data_t *out=NULL;
-  int ip=(flags & GAL_ARITHMETIC_FLAG_INPLACE) || (flags & 
GAL_ARITHMETIC_FLAG_FREE);
+  int ip= (    (flags & GAL_ARITHMETIC_FLAG_INPLACE)
+            || (flags & GAL_ARITHMETIC_FLAG_FREE)    );
+
+  /* The dataset may be empty. In this case, the output should also be
+     empty (we can have tables and images with 0 rows or pixels!). */
+  if(input->size==0 || input->array==NULL) return input;
 
+  /* Do the operation. */
   switch(operator)
     {
     case GAL_ARITHMETIC_OP_MINVAL:   out=gal_statistics_minimum(input);break;
@@ -673,7 +696,8 @@ arithmetic_from_statistics(int operator, int flags, 
gal_data_t *input)
 
 /* The size operator. Reports the size along a given dimension. */
 static gal_data_t *
-arithmetic_mknoise(int operator, int flags, gal_data_t *in, gal_data_t *arg)
+arithmetic_mknoise(int operator, int flags, gal_data_t *in,
+                   gal_data_t *arg)
 {
   gsl_rng *rng;
   const char *rng_name;
@@ -685,6 +709,10 @@ arithmetic_mknoise(int operator, int flags, gal_data_t 
*in, gal_data_t *arg)
      columns. */
   static unsigned long colcounter=0;
 
+  /* The dataset may be empty. In this case, the output should also be
+     empty (we can have tables and images with 0 rows or pixels!). */
+  if(in->size==0 || in->array==NULL) return in;
+
   /* Sanity checks. */
   if(arg->size!=1)
     error(EXIT_FAILURE, 0, "the first popped operand to the '%s' "
@@ -803,6 +831,7 @@ arithmetic_size(int operator, int flags, gal_data_t *in, 
gal_data_t *arg)
   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)
@@ -932,9 +961,23 @@ static void
 arithmetic_where(int flags, gal_data_t *out, gal_data_t *cond,
                  gal_data_t *iftrue)
 {
+  size_t i;
   int chb;    /* Read as: "Condition-Has-Blank" */
   unsigned char *c=cond->array, cb=GAL_BLANK_UINT8;
 
+  /* The datasets may be empty. In this case the output should also be
+     empty (we can have tables and images with 0 rows or pixels!). */
+  if(    out->size==0    || out->array==NULL
+      || cond->size==0   || cond->array==NULL
+      || iftrue->size==0 || iftrue->array==NULL )
+    {
+      if(flags & GAL_ARITHMETIC_FLAG_FREE)
+        { gal_data_free(cond); gal_data_free(iftrue); }
+      if(out->array) {free(out->array); out->array=NULL;}
+      if(out->dsize) for(i=0;i<out->ndim;++i) out->dsize[i]=0;
+      out->size=0; return;
+    }
+
   /* The condition operator has to be unsigned char. */
   if(cond->type!=GAL_TYPE_UINT8)
     error(EXIT_FAILURE, 0, "%s: the condition operand must be an "
@@ -1550,6 +1593,11 @@ arithmetic_multioperand(int operator, int flags, 
gal_data_t *list,
               "operator must be same", __func__,
               gal_arithmetic_operator_string(operator));
 
+      /* Make sure they actually have any data. */
+      if(tmp->size==0 || tmp->array==NULL)
+        error(EXIT_FAILURE, 0, "%s: atleast one input operand doesn't "
+              "have any data", __func__);
+
       /* Check the sizes. */
       if( gal_dimension_is_different(list, tmp) )
         error(EXIT_FAILURE, 0, "%s: the sizes of all operands to the '%s' "
@@ -1781,6 +1829,16 @@ arithmetic_binary(int operator, int flags, gal_data_t 
*l, gal_data_t *r)
   int quietmmap=l->quietmmap && r->quietmmap;
 
 
+  /* The datasets may be empty. In this case, the output should also be
+     empty (we can have tables and images with 0 rows or pixels!). */
+  if( l->size==0 || l->array==NULL || r->size==0 || r->array==NULL )
+    {
+      if(l->array==0 || l->array==NULL)
+        {   if(flags & GAL_ARITHMETIC_FLAG_FREE) gal_data_free(r); return l;}
+      else {if(flags & GAL_ARITHMETIC_FLAG_FREE) gal_data_free(l); return r;}
+    }
+
+
   /* Simple sanity check on the input sizes */
   if( !( (flags & GAL_ARITHMETIC_FLAG_NUMOK) && (l->size==1 || r->size==1))
       && gal_dimension_is_different(l, r) )
@@ -1948,6 +2006,17 @@ arithmetic_function_binary_flt(int operator, int flags, 
gal_data_t *il,
   double pi=3.14159265358979323846264338327;
   int quietmmap=il->quietmmap && ir->quietmmap;
 
+
+  /* The datasets may be empty. In this case, the output should also be
+     empty (we can have tables and images with 0 rows or pixels!). */
+  if( il->size==0 || il->array==NULL || ir->size==0 || ir->array==NULL )
+    {
+      if(il->array==0 || il->array==NULL)
+        {   if(flags & GAL_ARITHMETIC_FLAG_FREE) gal_data_free(ir); return il;}
+      else {if(flags & GAL_ARITHMETIC_FLAG_FREE) gal_data_free(il); return ir;}
+    }
+
+
   /* Simple sanity check on the input sizes */
   if( !( (flags & GAL_ARITHMETIC_FLAG_NUMOK) && (il->size==1 || ir->size==1))
       && gal_dimension_is_different(il, ir) )
@@ -2068,6 +2137,19 @@ arithmetic_box_around_ellipse(gal_data_t *d1, gal_data_t 
*d2,
           "function don't have the same number of elements",
           __func__);
 
+  /* The datasets may be empty. In this case the output should also be
+     empty (we can have tables and images with 0 rows or pixels!). */
+  if(    d1->size==0 || d1->array==NULL
+      || d2->size==0 || d2->array==NULL
+      || d3->size==0 || d3->array==NULL )
+    {
+      if(flags & GAL_ARITHMETIC_FLAG_FREE)
+        { gal_data_free(d2); gal_data_free(d3); }
+      if(d1->array) {free(d1->array); d1->array=NULL;}
+      if(d1->dsize) for(i=0;i<d1->ndim;++i) d1->dsize[i]=0;
+      d1->size=0; return d1;
+    }
+
   /* Convert the two inputs into double. Note that if the user doesn't want
      to free the inputs, we should make a copy of 'a_data' and 'b_data'
      because the output will also be written in them. */
@@ -2665,7 +2747,6 @@ gal_arithmetic(int operator, size_t numthreads, int 
flags, ...)
       d2 = va_arg(va, gal_data_t *);
       out=arithmetic_binary(operator, flags, d1, d2);
       break;
-
     case GAL_ARITHMETIC_OP_BITNOT:
       d1 = va_arg(va, gal_data_t *);
       out=arithmetic_bitwise_not(flags, d1);
diff --git a/lib/blank.c b/lib/blank.c
index 44c14fd8..187d4af3 100644
--- a/lib/blank.c
+++ b/lib/blank.c
@@ -563,6 +563,11 @@ gal_blank_flag(gal_data_t *input)
   gal_data_t *out;
   char **str=input->array, **strf=str+input->size;
 
+  /* The datasets may be empty. In this case the output should also be
+     empty (we can have tables and images with 0 rows or pixels!). */
+  if(input->size==0 || input->array==NULL) return input;
+
+  /* Do all the checks and allocations if a blank is actually present! */
   if( gal_blank_present(input, 0) )
     {
       /* Allocate a non-cleared output array, we are going to parse the
diff --git a/lib/data.c b/lib/data.c
index 2a2cbcec..ce4700c5 100644
--- a/lib/data.c
+++ b/lib/data.c
@@ -868,9 +868,22 @@ gal_data_copy_to_allocated(gal_data_t *in, gal_data_t *out)
   /* Correct the sizes of the output to be the same as the input. If it is
      equal, there is no problem, if not, the size information will be
      changed, so if you want to use this allocated space again, be sure to
-     re-set the size parameters. */
+     re-set the size parameters.
+
+     As defined in 'gal_data_initialize', when there is no dataset, 'dsize'
+     should be NULL. So if 'in->dsize==NULL', but 'out->dsize!=NULL', then
+     we should free 'out->dsize' and set it to NULL. Alternatively, if
+     'out' hasn't allocated a 'dsize', but 'in' already has one, allocate
+     'out->dsize', then fill it. */
   out->size=in->size;
-  memcpy(out->dsize, in->dsize, in->ndim * sizeof *(in->dsize) );
+  if(in->dsize)
+    {
+      if(out->dsize==NULL)
+        out->dsize=gal_pointer_allocate(GAL_TYPE_SIZE_T, in->ndim, 0,
+                                        __func__, "out->dsize");
+      memcpy(out->dsize, in->dsize, in->ndim * sizeof *(in->dsize) );
+    }
+  else if(out->dsize) { free(out->dsize); out->dsize=NULL; }
 }
 
 
diff --git a/lib/wcs.c b/lib/wcs.c
index 4fd0c90d..e29bbeea 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -2311,9 +2311,19 @@ gal_wcs_world_to_img(gal_data_t *coords, struct wcsprm 
*wcs, int inplace)
   int status, *stat=NULL, ncoord=coords->size, nelem;
   double *phi=NULL, *theta=NULL, *world=NULL, *pixcrd=NULL, *imgcrd=NULL;
 
+  /* It can happen that the input datasets are empty. In this case, simply
+     return them. */
+  if(coords->size==0 || coords->array==NULL)
+    {
+      if(inplace) return coords;
+      else error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at "
+                 "'%s' to fix the problem. The input has no data and "
+                 "'inplace' is not called", __func__, PACKAGE_BUGREPORT);
+    }
+
   /* Some sanity checks. */
-  wcs_convert_sanity_check_alloc(coords, wcs, __func__, &stat, &phi, &theta,
-                                 &world, &pixcrd, &imgcrd);
+  wcs_convert_sanity_check_alloc(coords, wcs, __func__, &stat, &phi,
+                                 &theta, &world, &pixcrd, &imgcrd);
   nelem=wcs->naxis; /* We have to make sure a WCS is given first. */
 
 
@@ -2323,7 +2333,8 @@ gal_wcs_world_to_img(gal_data_t *coords, struct wcsprm 
*wcs, int inplace)
 
 
   /* Use WCSLIB's wcss2p for the conversion. */
-  status=wcss2p(wcs, ncoord, nelem, world, phi, theta, imgcrd, pixcrd, stat);
+  status=wcss2p(wcs, ncoord, nelem, world, phi, theta, imgcrd,
+                pixcrd, stat);
   if(status)
     error(EXIT_FAILURE, 0, "%s: wcss2p ERROR %d: %s", __func__, status,
           wcs_errmsg[status]);