[gnuastro-commits] master 2f0d95a 4/4: Reading special floating point values in FITS ASCII tables

[Mohammad Akhlaghi]

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

    Reading special floating point values in FITS ASCII tables
    
    Until now, we would rely on CFITSIO to interpret the numerical values of
    FITS ASCII tables (as we do for the binary tables). However, there was a
    problem regarding special values like `INF' and `-INF'. In FITS ASCII
    tables, CFITSIO can't interpret these values (that it writes itself).
    
    To fix the problem, `fits_tab_read_ascii_float_special' was added. This
    function will read the column contents as strings, then it will parse them
    as floating point numbers with the robust `strtod' implementation of the
    GNU C Library (which is shipped with Gnuastro for internal usage).
    
    Also, the example code in `gal_fits_key_read_from_ptr' was corrected for a
    few mistakes.
---
 doc/gnuastro.texi |   9 ++--
 lib/fits.c        | 126 ++++++++++++++++++++++++++++++++++++++++++++----------
 2 files changed, 110 insertions(+), 25 deletions(-)

diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index ac4cc2c..fac52c5 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -22801,8 +22801,8 @@ Here is one example of using this function:
 /* Allocate an array of datasets. */
 gal_data_t *keysll=gal_data_array_calloc(N);
 
-/* Use the array as a list.*/
-for(i=0;i<N-2;++i) keysll[i].next=keysll[i+1];
+/* Make the array usable as a list too (by setting `next'). */
+for(i=0;i<N-1;++i) keysll[i].next=keysll[i+1];
 
 /* Fill the datasets with a `name' and a `type'. */
 keysll[0].name="NAME1";     keysll[0].type=GAL_TYPE_INT32;
@@ -22815,7 +22815,10 @@ gal_fits_key_read_from_ptr(fptr, keysll, 0, 0);
 
 /* Use the values as you like... */
 
-/* Free all the allocated spaces. */
+/* Free all the allocated spaces. Note that `name' wasn't allocated
+   in this example, so we should explicitly set it to NULL before
+   calling `gal_data_array_free'. */
+for(i=0;i<N;++i) keysll[i].name=NULL;
 gal_data_array_free(keysll, N, 1);
 @end example
 
diff --git a/lib/fits.c b/lib/fits.c
index 4e5021a..e584704 100644
--- a/lib/fits.c
+++ b/lib/fits.c
@@ -2343,6 +2343,65 @@ gal_fits_tab_info(char *filename, char *hdu, size_t 
*numcols,
 
 
 
+/* Read CFITSIO un-readable (INF, -INF or NAN) floating point values in
+   FITS ASCII tables. */
+static void
+fits_tab_read_ascii_float_special(char *filename, char *hdu, fitsfile *fptr,
+                                  gal_data_t *out, size_t colnum,
+                                  size_t numrows, size_t minmapsize)
+{
+  double tmp;
+  char **strarr;
+  char *tailptr;
+  gal_data_t *strrows;
+  size_t i, colwidth=50;
+  int anynul=0, status=0;
+
+  /* Allocate the dataset to keep the string values. */
+  strrows=gal_data_alloc(NULL, GAL_TYPE_STRING, 1, &numrows, NULL, 0,
+                         minmapsize, NULL, NULL, NULL);
+
+  /* Allocate the space to keep the string values. */
+  strarr=strrows->array;
+  for(i=0;i<numrows;++i)
+    {
+      errno=0;
+      strarr[i]=calloc(colwidth, sizeof *strarr[i]);
+      if(strarr[i]==NULL)
+        error(EXIT_FAILURE, errno, "%s: allocating %zu bytes for "
+              "strarr[%zu]", __func__, colwidth * sizeof *strarr[i], i);
+    }
+
+  /* Read the column as a string. */
+  fits_read_col(fptr, TSTRING, colnum, 1, 1, out->size, NULL, strrows->array,
+                &anynul, &status);
+  gal_fits_io_error(status, NULL);
+
+  /* Convert the strings to float. */
+  for(i=0;i<numrows;++i)
+    {
+      /* Parse the string. */
+      tmp=strtod(strarr[i], &tailptr);
+      if(tailptr==strarr[i])
+        error(EXIT_FAILURE, 0, "%s (hdu %s): couldn't parse row %zu of "
+              "column %zu (with value `%s') as a floating point number",
+              filename, hdu, i+1, colnum, strarr[i]);
+
+      /* Write it into the output dataset. */
+      if(out->type==GAL_TYPE_FLOAT32)
+        ((float *)(out->array))[i]=tmp;
+      else
+        ((double *)(out->array))[i]=tmp;
+    }
+
+  /* Clean up. */
+  gal_data_free(strrows);
+}
+
+
+
+
+
 /* Read the column indexs given in the `indexll' linked list from a FITS
    table into a linked list of data structures, note that this is a
    low-level function, so the output data linked list is the inverse of the
@@ -2358,12 +2417,17 @@ gal_fits_tab_read(char *filename, char *hdu, size_t 
numrows,
   char **strarr;
   fitsfile *fptr;
   gal_data_t *out=NULL;
-  int status=0, anynul=0;
   gal_list_sizet_t *ind;
+  int isfloat, status=0, anynul=0, hdutype;
 
   /* Open the FITS file */
   fptr=gal_fits_hdu_open_format(filename, hdu, 1);
 
+  /* See if its a Binary or ASCII table (necessary for floating point blank
+     values). */
+  if (fits_get_hdu_type(fptr, &hdutype, &status) )
+    gal_fits_io_error(status, NULL);
+
   /* Pop each index and read/store the array. */
   for(ind=indexll; ind!=NULL; ind=ind->next)
     {
@@ -2379,35 +2443,53 @@ gal_fits_tab_read(char *filename, char *hdu, size_t 
numrows,
          disp_width element of the data structure, which is done
          automatically in `gal_fits_table_info'. */
       if(out->type==GAL_TYPE_STRING)
-        for(i=0;i<numrows;++i)
-          {
-            strarr=out->array;
-            errno=0;
-            strarr[i]=calloc(allcols[ind->v].disp_width+1, sizeof *strarr[i]);
-            if(strarr[i]==NULL)
-              error(EXIT_FAILURE, errno, "%s: allocating %zu bytes for "
-                    "strarr[%zu]", __func__,
-                    (allcols[ind->v].disp_width+1) * sizeof *strarr[i], i);
-          }
+        {
+          strarr=out->array;
+          for(i=0;i<numrows;++i)
+            {
+              errno=0;
+              strarr[i]=calloc(allcols[ind->v].disp_width+1,
+                               sizeof *strarr[i]);
+              if(strarr[i]==NULL)
+                error(EXIT_FAILURE, errno, "%s: allocating %zu bytes for "
+                      "strarr[%zu]", __func__,
+                      (allcols[ind->v].disp_width+1) * sizeof *strarr[i], i);
+            }
+        }
 
       /* Allocate a blank value for the given type and read/store the
-         column using CFITSIO. Afterwards, free the blank value. Note that
-         we only need blank values for integer types. For floating point
-         types, the FITS standard defines blanks as NaN (same as almost any
-         other software like Gnuastro). However if a blank value is
-         specified, CFITSIO will convert other special numbers like `inf'
-         to NaN also. We want to be able to distringuish `inf' and NaN
-         here, so for floating point types, we won't define any blank
-         value. */
-      blank = ( (out->type==GAL_TYPE_FLOAT32 || out->type==GAL_TYPE_FLOAT64)
+         column using CFITSIO. Note that for binary tables, we only need
+         blank values for integer types. For binary floating point types,
+         the FITS standard defines blanks as NaN (same as almost any other
+         software like Gnuastro). However if a blank value is specified,
+         CFITSIO will convert other special numbers like `inf' to NaN
+         also. We want to be able to distringuish `inf' and NaN here, so
+         for floating point types in binary tables, we won't define any
+         blank value. In ASCII tables, CFITSIO doesn't read the `NAN'
+         values (that it has written itself) unless we specify a blank
+         pointer/value. */
+      isfloat = out->type==GAL_TYPE_FLOAT32 || out->type==GAL_TYPE_FLOAT64;
+      blank = ( ( hdutype==BINARY_TBL && isfloat )
                 ? NULL
                 : gal_blank_alloc_write(out->type) );
       fits_read_col(fptr, gal_fits_type_to_datatype(out->type), ind->v+1,
                     1, 1, out->size, blank, out->array, &anynul, &status);
-      gal_fits_io_error(status, NULL);
 
-      /* Clean up. */
+      /* In the ASCII table format, CFITSIO might not be able to read `INF'
+         or `-INF'. In this case, it will set status to `BAD_C2D' or
+         `BAD_C2F'. So, we'll use our own parser for the column values. */
+      if( hdutype==ASCII_TBL
+          && isfloat
+          && (status==BAD_C2D || status==BAD_C2F) )
+        {
+          fits_tab_read_ascii_float_special(filename, hdu, fptr, out,
+                                            ind->v+1, numrows, minmapsize);
+          status=0;
+        }
+
+      /* Clean up and sanity check. */
       if(blank) free(blank);
+      gal_fits_io_error(status, NULL);
     }
 
   /* Close the FITS file */