[gnuastro-commits] master 7356f18 2/2: Fixed 32-bit check errors and freeing empty WCS

[Mohammad Akhlaghi]

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

    Fixed 32-bit check errors and freeing empty WCS
    
    After Gnuastro 0.3.7 was uploaded to the Debian server and tested on many
    systems, it was found that two tests fail (one in CosmicCalculator) and
    another in MakeCatalog. So after installing an i386 Debian as a virtual
    machine and testing it, the causes of the problem were found and corrected:
    
     - The CosmicCalculator test was corrected by defining a `sum' variable to
       be the sum of the different densities and using that for a check. I
       guess this is a compiler issue. But using that variable, also helped in
       several places where this sum was necessary in the function. Also, since
       users might use scripts to generate the densities, we now allow for
       small floating point errors also.
    
     - The problem in MakeCatalog was due to not being able to read the
       `NUMLABS' keyword in the clumps image. The problem was that on 32-bit
       systems, `long' and `int' have the same length of 32-bits. But the
       CFITSIO manual is silent on the length of its `TINT' (or `TUINT')
       lengths. So on a 32-bit system, when trying to read a number into
       `TUINT' it would give an overflow error. As a result, when converting
       Gnuastro's types to CFITSIO's types, for 32-bit integers, CFITSIO's LONG
       types are preferred.
    
    Also with this commit, the work done in `gal_wcs_read_fitsptr' of the last
    commit was cleaned and some bugs in it were corrected, for example the
    `altlin' variable is now only checked when a WCS structure was actually
    read.
---
 bin/cosmiccal/ui.c     | 16 ++++++++-------
 bin/mkcatalog/ui.c     |  8 ++++++--
 lib/fits.c             | 17 +++++++++------
 lib/wcs.c              | 56 +++++++++++++++++++++++++++++++++++---------------
 tests/arithmetic/or.sh |  2 +-
 5 files changed, 67 insertions(+), 32 deletions(-)

diff --git a/bin/cosmiccal/ui.c b/bin/cosmiccal/ui.c
index c2864de..a7d8d98 100644
--- a/bin/cosmiccal/ui.c
+++ b/bin/cosmiccal/ui.c
@@ -195,16 +195,18 @@ parse_opt(int key, char *arg, struct argp_state *state)
 static void
 ui_read_check_only_options(struct cosmiccalparams *p)
 {
-  /* Check if the density fractions add up to 1. */
-  if( (p->olambda + p->omatter + p->oradiation) != 1.0f )
-    error(EXIT_FAILURE, 0, "sum of fractional densities is not 1, but %f. "
+  double sum=p->olambda + p->omatter + p->oradiation;
+
+  /* Check if the density fractions add up to 1 (within floating point
+     error). */
+  if( sum > (1+1e-8) || sum < (1-1e-8) )
+    error(EXIT_FAILURE, 0, "sum of fractional densities is not 1, but %.8f. "
           "The cosmological constant (`olambda'), matter (`omatter') "
-          "and radiation (`oradiation') densities are given as %f, %f, %f",
-          p->olambda + p->omatter + p->oradiation, p->olambda, p->omatter,
-          p->oradiation);
+          "and radiation (`oradiation') densities are given as %.8f, %.8f, "
+          "%.8f", sum, p->olambda, p->omatter, p->oradiation);
 
   /* The curvature fractional density: */
-  p->ocurv=1-(p->olambda+p->omatter+p->oradiation);
+  p->ocurv=1-sum;
 
   /* Convert H0 from km/sec/Mpc to 1/sec: */
   p->H0s=p->H0/1000/GSL_CONST_MKSA_PARSEC;
diff --git a/bin/mkcatalog/ui.c b/bin/mkcatalog/ui.c
index 7096feb..87f95a1 100644
--- a/bin/mkcatalog/ui.c
+++ b/bin/mkcatalog/ui.c
@@ -500,6 +500,7 @@ ui_preparations_read_inputs(struct mkcatalogparams *p)
 static void
 ui_preparations_read_keywords(struct mkcatalogparams *p)
 {
+  char *msg;
   gal_data_t *tmp;
   gal_data_t *keys=gal_data_array_calloc(2);
   char *stdfile=p->stdfile ? p->stdfile : p->inputname;
@@ -575,8 +576,11 @@ ui_preparations_read_keywords(struct mkcatalogparams *p)
       gal_fits_key_read(clumpsfile, p->clumpshdu, keys, 0, 0);
       if(keys[0].status) p->clumpsn=NAN;
       if(keys[1].status)
-        error(EXIT_FAILURE, 0, "couldn't find NCLUMPS in the header of "
-              "%s (hdu: %s).", p->clumpsfile, p->clumpshdu);
+        {
+          asprintf(&msg, "couldn't find/read NUMLABS in the header of "
+                   "%s (hdu: %s), see error above", clumpsfile, p->clumpshdu);
+          gal_fits_io_error(keys[1].status, msg);
+        }
     }
 
 
diff --git a/lib/fits.c b/lib/fits.c
index d5f555d..ef1f328 100644
--- a/lib/fits.c
+++ b/lib/fits.c
@@ -308,17 +308,22 @@ gal_fits_type_to_datatype(uint8_t type)
       else if( sizeof(int)      == w )   return TINT;
       break;
 
+    /* On 32-bit systems, the length of `int' and `long' are both
+       32-bits. But CFITSIO's LONG type is preferred because it is designed
+       to be 32-bit. Its `INT' type is not clearly defined and caused
+       problems when reading keywords.*/
     case GAL_TYPE_UINT32:
       w=4;
-      if     ( sizeof(int)      == w )   return TUINT;
-      else if( sizeof(long)     == w )   return TULONG;
+      if     ( sizeof(long)     == w )   return TULONG;
+      else if( sizeof(int)      == w )   return TUINT;
       else if( sizeof(short)    == w )   return TUSHORT;
       break;
 
+    /* Similar to UINT32 above. */
     case GAL_TYPE_INT32:
       w=4;
-      if     ( sizeof(int)      == w )   return TINT;
-      else if( sizeof(long)     == w )   return TLONG;
+      if     ( sizeof(long)     == w )   return TLONG;
+      else if( sizeof(int)      == w )   return TINT;
       else if( sizeof(short)    == w )   return TSHORT;
       break;
 
@@ -1322,8 +1327,8 @@ gal_fits_img_info(fitsfile *fptr, int *type, size_t 
*ndim, size_t **dsize,
         {
         switch(i)
           {
-          case 4: if(unit) {str = key->array; *unit = *str;}   break;
-          case 3: if(name) {str = key->array; *name = *str;}   break;
+          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;
           default:
diff --git a/lib/wcs.c b/lib/wcs.c
index eb566e9..f99eea0 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -70,8 +70,8 @@ gal_wcs_read_fitsptr(fitsfile *fptr, size_t hstartwcs, size_t 
hendwcs,
                      int *nwcs)
 {
   /* Declaratins: */
-  struct wcsprm *wcs;
   int nkeys=0, status=0;
+  struct wcsprm *wcs=NULL;
   char *fullheader, *to, *from;
   int relax    = WCSHDR_all; /* Macro: use all informal WCS extensions. */
   int ctrl     = 0;          /* Don't report why a keyword wasn't used. */
@@ -111,7 +111,8 @@ gal_wcs_read_fitsptr(fitsfile *fptr, size_t hstartwcs, 
size_t hendwcs,
       /*******************************************************/
     }
 
-  /* WCSlib function */
+
+  /* WCSlib function to parse the FITS headers. */
   status=wcspih(fullheader, nkeys, relax, ctrl, &nreject, nwcs, &wcs);
   if(status)
     {
@@ -125,22 +126,45 @@ gal_wcs_read_fitsptr(fitsfile *fptr, size_t hstartwcs, 
size_t hendwcs,
     gal_fits_io_error(status, "problem in fitsarrayvv.c for freeing "
                            "the memory used to keep all the headers");
 
+
   /* Set the internal structure: */
-  status=wcsset(wcs);
-  if(status)
+  if(wcs)
     {
-      fprintf(stderr, "\n##################\n"
-              "WCSLIB Warning: wcsset ERROR %d: %s.\n"
-              "##################\n",
-            status, wcs_errmsg[status]);
-      wcs=NULL; *nwcs=0;
+      /* CTYPE is a mandatory WCS keyword, so if it hasn't been given (its
+         '\0'), then the headers didn't have a WCS structure. However,
+         WCSLIB still fills in the basic information (for example the
+         dimensionality of the dataset). */
+      if(wcs->ctype[0][0]=='\0')
+        {
+          wcsfree(wcs);
+          wcs=NULL;
+          *nwcs=0;
+        }
+      else
+        {
+          /* Set the WCS structure. */
+          status=wcsset(wcs);
+          if(status)
+            {
+              fprintf(stderr, "\n##################\n"
+                      "WCSLIB Warning: wcsset ERROR %d: %s.\n"
+                      "##################\n",
+                      status, wcs_errmsg[status]);
+              wcsfree(wcs);
+              wcs=NULL;
+              *nwcs=0;
+            }
+          else
+            /* A correctly useful WCS is present. When no PC matrix
+               elements were present in the header, the default PC matrix
+               (a unity matrix) is used. In this case WCSLIB doesn't set
+               `altlin' (and gives it a value of 0). In Gnuastro, later on,
+               we might need to know the type of the matrix used, so in
+               such a case, we will set `altlin' to 1. */
+            if(wcs->altlin==0) wcs->altlin=1;
+        }
     }
 
-  /* Initialize the wcsprm struct
-  if ((status = wcsset(*wcs)))
-    error(EXIT_FAILURE, 0, "%s: wcsset ERROR %d: %s.\n", __func__,
-          status, wcs_errmsg[status]);
-  */
 
   /* Return the WCS structure. */
   return wcs;
@@ -294,13 +318,13 @@ gal_wcs_warp_matrix(struct wcsprm *wcs)
           __func__, size*sizeof *out);
 
   /* Fill in the array. */
-  if(wcs->altlin & 0x1)        /* Has a PCi_j array. */
+  if(wcs->altlin & 0x1)          /* Has a PCi_j array. */
     {
       for(i=0;i<wcs->naxis;++i)
         for(j=0;j<wcs->naxis;++j)
           out[i*wcs->naxis+j] = wcs->cdelt[i] * wcs->pc[i*wcs->naxis+j];
     }
-  else if(wcs->altlin & 0x2)     /* Has CDi_j array */
+  else if(wcs->altlin & 0x2)     /* Has CDi_j array.   */
     {
       for(i=0;i<size;++i)
         out[i]=wcs->cd[i];
diff --git a/tests/arithmetic/or.sh b/tests/arithmetic/or.sh
index eb08426..263dcbf 100755
--- a/tests/arithmetic/or.sh
+++ b/tests/arithmetic/or.sh
@@ -49,4 +49,4 @@ if [ ! -f $img      ]; then echo "$img does not exist.";   
exit 77; fi
 
 # Actual test script
 # ==================
-$execname $img 6 eq $img 3 eq or -h1 -h1 --output=or.fits
+$execname $img 6 eq $img 3 eq or -h2 -h2 --output=or.fits