[gnuastro-commits] master 6b53c5e 015/113: Bug fixes in master: commits 1ff1c25 and 540af65

[Mohammad Akhlaghi]

branch: master
commit 6b53c5ee0a234afbd70cccf6726c202f0d75d41d
Author: Mohammad Akhlaghi <akhlaghi@gnu.org>
Commit: Mohammad Akhlaghi <akhlaghi@gnu.org>

    Bug fixes in master: commits 1ff1c25 and 540af65
    
    Two recent bugs were found in Gnuastro and fixed in the master branch. With
    this commit, they are also copied/implemented here. I didn't actually use
    Git's merge because the history can get very messy with all the merges with
    master. It is cleaner to manually copy the files like this.
---
 NEWS                        |  3 ++
 bin/arithmetic/arithmetic.c | 28 ++++++++++++-
 bin/arithmetic/operands.c   |  6 +--
 lib/fits.c                  | 99 ++++++++++++++++++++++++++++++---------------
 4 files changed, 98 insertions(+), 38 deletions(-)

diff --git a/NEWS b/NEWS
index 64cff8e..c59f194 100644
--- a/NEWS
+++ b/NEWS
@@ -99,6 +99,9 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
 
   Crashes on 32-bit and big-endian systems (bug #51476).
 
+  Reading BZERO for unsigned 64-bit integers (bug #51555).
+
+  Arithmetic with one file and no operators (bug #51559).
 
 
 
diff --git a/bin/arithmetic/arithmetic.c b/bin/arithmetic/arithmetic.c
index 1456a12..a82e256 100644
--- a/bin/arithmetic/arithmetic.c
+++ b/bin/arithmetic/arithmetic.c
@@ -29,6 +29,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <string.h>
 #include <stdlib.h>
 
+#include <gnuastro/wcs.h>
 #include <gnuastro/fits.h>
 #include <gnuastro/arithmetic.h>
 
@@ -132,6 +133,7 @@ void
 reversepolish(struct arithmeticparams *p)
 {
   int op=0, nop=0;
+  char *filename, *hdu;
   unsigned int numop, i;
   gal_list_str_t *token;
   gal_data_t *d1=NULL, *d2=NULL, *d3=NULL;
@@ -333,12 +335,34 @@ reversepolish(struct arithmeticparams *p)
     error(EXIT_FAILURE, 0, "too many operands");
 
 
-  /* Set the output type. */
-  d1=p->operands->data;
+  /* If the final operand has a filename, but it `data' element is NULL,
+     then the file hasn't actually be read yet. In this case, we need to
+     read the contents of the file and put the resulting dataset into the
+     operands `data' element. This can happen for example if no operators
+     are called and there is only one filename as an argument (which can
+     happen in scripts). */
+  if(p->operands->data==NULL && p->operands->filename)
+    {
+      /* Read the desired image and report it if necessary. */
+      hdu=p->operands->hdu;
+      filename=p->operands->filename;
+      if( gal_fits_name_is_fits(filename) )
+        {
+          p->operands->data=gal_fits_img_read(filename,hdu,p->cp.minmapsize);
+          p->refdata.wcs=gal_wcs_read(filename, hdu, 0, 0, &p->refdata.nwcs);
+          if(!p->cp.quiet) printf(" - %s (hdu %s) is read.\n", filename, hdu);
+        }
+      else
+        error(EXIT_FAILURE, 0, "%s: a bug! please contact us at %s to fix "
+              "the problem. While `operands->data' is NULL, the filename "
+              "(`%s') is not recognized as a FITS file", __func__,
+              PACKAGE_BUGREPORT, filename);
+    }
 
 
   /* If the final data structure has more than one element, write it as a
      FITS file. Otherwise, print it in the standard output. */
+  d1=p->operands->data;
   if(d1->size==1)
     {
       /* To simplify the printing process, we will first change it to
diff --git a/bin/arithmetic/operands.c b/bin/arithmetic/operands.c
index c069f34..585ea0b 100644
--- a/bin/arithmetic/operands.c
+++ b/bin/arithmetic/operands.c
@@ -110,8 +110,8 @@ operands_pop(struct arithmeticparams *p, char *operator)
   /* If the operand linked list has finished, then give an error and
      exit. */
   if(operands==NULL)
-    error(EXIT_FAILURE, 0, "not enough operands for the \"%s\" "
-          "operator", operator);
+    error(EXIT_FAILURE, 0, "not enough operands for the \"%s\" operator",
+          operator);
 
 
   /* Set the dataset. If filename is present then read the file
@@ -126,7 +126,7 @@ operands_pop(struct arithmeticparams *p, char *operator)
       if(p->popcounter==0)
         p->refdata.wcs=gal_wcs_read(filename, hdu, 0, 0, &p->refdata.nwcs);
 
-      /* Read the input image as a double type array. */
+      /* Read the dataset. */
       data=gal_fits_img_read(filename, hdu, p->cp.minmapsize);
 
       /* When the reference data structure's dimensionality is non-zero, it
diff --git a/lib/fits.c b/lib/fits.c
index 2f3860d..6cf7621 100644
--- a/lib/fits.c
+++ b/lib/fits.c
@@ -471,39 +471,72 @@ gal_fits_datatype_to_type(int datatype, int 
is_table_column)
    different from the type from BITPIX. This function does the necessary
    corrections.*/
 static void
-fits_type_correct(int *type, double bscale, double bzero)
+fits_type_correct(int *type, double bscale, char *bzero_str)
 {
+  double bzero;
   int tofloat=1;
+  char *tailptr, *bzero_u64="9223372036854775808";
 
-  /* Work based on type. For the default conversions defined by the FITS
-     standard to change the signs of integers, make the proper correction,
-     otherwise set the type to float. */
-  if(bscale==1.0f)
-    switch(*type)
-      {
-      case GAL_TYPE_UINT8:
-        if(bzero == -128)         { *type = GAL_TYPE_INT8;   tofloat=0; }
-        break;
+  /* If bzero_str is given and `bscale=1.0' the case might be that we are
+     dealing with an integer dataset that just needs a different sign. */
+  if(bzero_str && bscale==1.0f)
+    {
+      /* Read the `bzero' string as a `double' number. */
+      bzero  = strtod(bzero_str, &tailptr);
+      if(tailptr==bzero_str)
+        error(EXIT_FAILURE, 0, "%s: BZERO value `%s' couldn't be "
+              "parsed as a number", __func__, bzero_str);
+
+      /* Work based on type. For the default conversions defined by the
+         FITS standard to change the signs of integers, make the proper
+         correction, otherwise set the type to float. */
+      switch(*type)
+        {
+        case GAL_TYPE_UINT8:
+          if(bzero == -128.0f)      { *type = GAL_TYPE_INT8;   tofloat=0; }
+          break;
 
-      case GAL_TYPE_INT16:
-        if(bzero == 32768)        { *type = GAL_TYPE_UINT16; tofloat=0; }
-        break;
+        case GAL_TYPE_INT16:
+          if(bzero == 32768)        { *type = GAL_TYPE_UINT16; tofloat=0; }
+          break;
 
-      case GAL_TYPE_INT32:
-        if(bzero == 2147483648LU) { *type = GAL_TYPE_UINT32; tofloat=0; }
-        break;
+        case GAL_TYPE_INT32:
+          if(bzero == 2147483648LU) { *type = GAL_TYPE_UINT32; tofloat=0; }
+          break;
 
-      case GAL_TYPE_INT64:
-        if(bzero == 9223372036854775808LLU)
-          {*type = GAL_TYPE_UINT64; tofloat=0;}
-        break;
+        /* The `bzero' value for unsigned 64-bit integers has 19 decimal
+           digits, but a 64-bit floating point (`double' type) can only
+           safely store 15 decimal digits. As a result, the safest way to
+           check the `bzero' value for this type is to compare it as a
+           string. But all integers nearby (for example
+           `9223372036854775807') get rounded to this same value (when
+           stored as `double'). So we will also check the parsed number and
+           if it equals this number, a warning will be printed. */
+        case GAL_TYPE_INT64:
+          if( !strcmp(bzero_str, bzero_u64) )
+            {*type = GAL_TYPE_UINT64; tofloat=0;}
+          else
+            if( bzero == 9223372036854775808LLU )
+              {
+                fprintf(stderr, "\nWARNING in %s: the BZERO header keyword "
+                        "value (`%s') is very close (but not exactly equal) "
+                        "to `%s'. The latter value in the FITS standard is "
+                        "used to identify that the dataset should be read as "
+                        "unsigned 64-bit integers instead of signed 64-bit "
+                        "integers. Depending on your version of CFITSIO, "
+                        "it may be read as a signed or unsigned 64-bit "
+                        "integer array\n\n", __func__, bzero_str, bzero_u64);
+                tofloat=0;
+              }
+          break;
 
-        /* For the other types (when `BSCALE=1.0f'), currently no correction is
-           necessary, maybe later we can check if the scales are integers and
-           set the integer output type to the smallest type that can allow the
-           scaled values. */
-      default: tofloat=0;
-      }
+          /* For the other types (when `BSCALE=1.0f'), currently no
+             correction is necessary, maybe later we can check if the
+             scales are integers and set the integer output type to the
+             smallest type that can allow the scaled values. */
+        default: tofloat=0;
+        }
+    }
 
   /* If the type must be a float, then do the conversion. */
   if(tofloat) *type=GAL_TYPE_FLOAT32;
@@ -1288,10 +1321,10 @@ void
 gal_fits_img_info(fitsfile *fptr, int *type, size_t *ndim, size_t **dsize,
                   char **name, char **unit)
 {
-  char **str;
+  double bscale=NAN;
   size_t i, dsize_key=1;
+  char **str, *bzero_str=NULL;
   int bitpix, status=0, naxis;
-  double bzero=NAN, bscale=NAN;
   gal_data_t *key, *keysll=NULL;
   long naxes[GAL_FITS_MAX_NDIM];
 
@@ -1314,7 +1347,7 @@ gal_fits_img_info(fitsfile *fptr, int *type, size_t 
*ndim, size_t **dsize,
                           NULL, 0, -1, "EXTNAME", NULL, NULL);
   gal_list_data_add_alloc(&keysll, NULL, GAL_TYPE_FLOAT64, 1, &dsize_key,
                           NULL, 0, -1, "BSCALE", NULL, NULL);
-  gal_list_data_add_alloc(&keysll, NULL, GAL_TYPE_FLOAT64, 1, &dsize_key,
+  gal_list_data_add_alloc(&keysll, NULL, GAL_TYPE_STRING, 1, &dsize_key,
                           NULL, 0, -1, "BZERO", NULL, NULL);
   gal_fits_key_read_from_ptr(fptr, keysll, 0, 0);
 
@@ -1331,8 +1364,8 @@ gal_fits_img_info(fitsfile *fptr, int *type, size_t 
*ndim, size_t **dsize,
           {
           case 4: if(unit) {str = key->array; *unit = *str; *str=NULL;} break;
           case 3: if(name) {str = key->array; *name = *str; *str=NULL;} break;
-          case 2: bscale = *(double *)(key->array);    break;
-          case 1: bzero  = *(double *)(key->array);    break;
+          case 2: bscale = *(double *)(key->array);                     break;
+          case 1: str = key->array; bzero_str = *str;                   break;
           default:
             error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
                   "fix the problem. For some reason, there are more "
@@ -1342,8 +1375,8 @@ gal_fits_img_info(fitsfile *fptr, int *type, size_t 
*ndim, size_t **dsize,
       ++i;
     }
 
-  if( !isnan(bscale) || !isnan(bzero) )
-    fits_type_correct(type, bscale, bzero);
+  if( !isnan(bscale) || bzero_str )
+    fits_type_correct(type, bscale, bzero_str);
 
 
   /* Allocate the array to keep the dimension size and fill it in, note