[gnuastro-commits] master 2071999 08/16: Blank CFITSIO pointers set by datatype, not bitpix

[Mohammad Akhlaghi]

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

    Blank CFITSIO pointers set by datatype, not bitpix
    
    While working on the new Table utility, it was necessary to read a larger
    collection of datatypes than those of the BITPIX keyword of FITS
    images. CFITSIO defines its internal datatypes (starting with a `T') for
    this purpose. Since the table columns and images have some types in common,
    the more low-level parameter that should be used to set blank values and
    allocate space should be the datatypes, not image or table type keywords
    (which have different macros).
    
    In this commit, the `gal_fits_bitpix_blank' function was replaced with
    `gal_fits_datatype_blank' and now defines blank values for all CFITSIO
    datatypes in its manual. It uses the GNU Scientific Library (GSL) complex
    types for the complex datatypes.
---
 lib/fits.c          |  127 +++++++++++++++++++++++++++++++++++++++++----------
 lib/gnuastro/fits.h |   32 ++++++++++---
 src/imgcrop/crop.c  |   32 ++++++-------
 src/imgcrop/ui.c    |    2 +-
 4 files changed, 145 insertions(+), 48 deletions(-)

diff --git a/lib/fits.c b/lib/fits.c
index c4091a5..4b4dc6b 100644
--- a/lib/fits.c
+++ b/lib/fits.c
@@ -179,63 +179,142 @@ gal_fits_bitpix_to_dtype(int bitpix)
 
 
 void *
-gal_fits_bitpix_blank(int bitpix)
+gal_fits_datatype_blank(int datatype)
 {
+  /* Define the pointers, note that we are ordering them based on the
+     CFITSIO manual to be more easily comparable. */
   unsigned char *b;
+  char *c;
+  char **str;
   short *s;
   long *l;
   LONGLONG *L;
   float *f;
   double *d;
+  gsl_complex_float *complex;
+  gsl_complex *dblcomplex;
+  int *i;
+  unsigned int *ui;
+  unsigned short *us;
 
   errno=0;
-  switch(bitpix)
+  switch(datatype)
     {
-    case BYTE_IMG:
-      b=malloc(sizeof(unsigned char));
+    case TBIT:
+      error(EXIT_FAILURE, 0, "Currently GSL doesn't support TBIT datatype, "
+            "please get in touch with us to implement it.");
+
+    case TBYTE:
+      b=malloc(sizeof *b);
       if(b==NULL)
-        error(EXIT_FAILURE, errno, "%lu bytes", sizeof(unsigned char));
+        error(EXIT_FAILURE, errno, "%lu bytes for blank TBYTE",
+              sizeof *b);
       *b=GAL_FITS_BYTE_BLANK;
       return b;
 
-    case SHORT_IMG:
-      s=malloc(sizeof(short));
+      /* CFITSIO says "int for keywords, char for table columns". Here we
+         are only assuming table columns. So in practice this also applies
+         to TSBYTE.*/
+    case TLOGICAL: case TSBYTE:
+      c=malloc(sizeof *c);
+      if(c==NULL)
+        error(EXIT_FAILURE, errno, "%lu bytes for blank TLOGICAL, or TSBYTE",
+              sizeof *c);
+      *c=GAL_FITS_LOGICAL_BLANK;
+      return c;
+
+    case TSTRING:
+      str=malloc(sizeof *str);
+      if(str==NULL)
+        error(EXIT_FAILURE, errno, "%lu bytes for blank TSTRING",
+              sizeof *s);
+      *str=GAL_FITS_STRING_BLANK;
+      return str;
+
+    case TSHORT:
+      s=malloc(sizeof *s);
       if(s==NULL)
-        error(EXIT_FAILURE, errno, "%lu bytes", sizeof(short));
+        error(EXIT_FAILURE, errno, "%lu bytes for blank TSHORT",
+              sizeof *s);
       *s=GAL_FITS_SHORT_BLANK;
       return s;
 
-    case LONG_IMG:
-      l=malloc(sizeof(long));
+    case TLONG:
+      l=malloc(sizeof *l);
       if(l==NULL)
-        error(EXIT_FAILURE, errno, "%lu bytes", sizeof(long));
+        error(EXIT_FAILURE, errno, "%lu bytes for blank TLONG",
+              sizeof *l);
       *l=GAL_FITS_LONG_BLANK;
       return l;
 
-    case LONGLONG_IMG:
-      L=malloc(sizeof(LONGLONG));
+    case TLONGLONG:
+      L=malloc(sizeof *L);
       if(L==NULL)
-        error(EXIT_FAILURE, errno, "%lu bytes", sizeof(LONGLONG));
+        error(EXIT_FAILURE, errno, "%lu bytes for blank TLONGLONG",
+              sizeof *L);
       *L=GAL_FITS_LLONG_BLANK;
       return L;
 
-    case FLOAT_IMG:
-      f=malloc(sizeof(float));
+    case TFLOAT:
+      f=malloc(sizeof *f);
       if(f==NULL)
-        error(EXIT_FAILURE, errno, "%lu bytes", sizeof(float));
+        error(EXIT_FAILURE, errno, "%lu bytes for blank TFLOAT",
+              sizeof *f);
       *f=GAL_FITS_FLOAT_BLANK;
       return f;
 
-    case DOUBLE_IMG:
-      d=malloc(sizeof(double));
+    case TDOUBLE:
+      d=malloc(sizeof *d);
       if(d==NULL)
-        error(EXIT_FAILURE, errno, "%lu bytes", sizeof(double));
+        error(EXIT_FAILURE, errno, "%lu bytes for blank TDOUBLE",
+              sizeof *d);
       *d=GAL_FITS_FLOAT_BLANK;
       return d;
 
+    case TCOMPLEX:
+      complex=malloc(sizeof *complex);
+      if(complex==NULL)
+        error(EXIT_FAILURE, errno, "%lu bytes for blank TCOMPLEX",
+              sizeof *complex);
+      GAL_FITS_TCOMPLEX_BLANK(complex);
+      return complex;
+
+    case TDBLCOMPLEX:
+      dblcomplex=malloc(sizeof *dblcomplex);
+      if(dblcomplex==NULL)
+        error(EXIT_FAILURE, errno, "%lu bytes for blank TDBLCOMPLEX",
+              sizeof *dblcomplex);
+      GAL_FITS_TCOMPLEX_BLANK(dblcomplex);
+      return dblcomplex;
+
+    case TINT:
+      i=malloc(sizeof *i);
+      if(i==NULL)
+        error(EXIT_FAILURE, errno, "%lu bytes for blank TINT",
+              sizeof *i);
+      *i=GAL_FITS_INT_BLANK;
+      return i;
+
+    case TUINT:
+      ui=malloc(sizeof *ui);
+      if(ui==NULL)
+        error(EXIT_FAILURE, errno, "%lu bytes for blank TUINT",
+              sizeof *ui);
+      *ui=GAL_FITS_UINT_BLANK;
+      return ui;
+
+    case TUSHORT:
+      us=malloc(sizeof *us);
+      if(us==NULL)
+        error(EXIT_FAILURE, errno, "%lu bytes for blank TUSHORT",
+              sizeof *us);
+      *ui=GAL_FITS_USHORT_BLANK;
+      return us;
+
+
     default:
-      error(EXIT_FAILURE, 0, "bitpix value of %d not recognized",
-            bitpix);
+      error(EXIT_FAILURE, 0, "datatype value of %d not recognized",
+            datatype);
     }
 
   return NULL;
@@ -1361,7 +1440,7 @@ gal_fits_hdu_to_array(char *filename, char *hdu, int 
*bitpix,
   *s1=naxes[0];
 
   /* Allocate space for the array. */
-  bitblank=gal_fits_bitpix_blank(*bitpix);
+  bitblank=gal_fits_datatype_blank( gal_fits_bitpix_to_dtype(*bitpix) );
   *array=gal_fits_bitpix_alloc(*s0 * *s1, *bitpix);
 
   /* Read the image into the allocated array: */
@@ -1429,7 +1508,7 @@ gal_fits_array_to_file(char *filename, char *hdu, int 
bitpix,
     if(bitpix==BYTE_IMG || bitpix==SHORT_IMG
        || bitpix==LONG_IMG || bitpix==LONGLONG_IMG)
       {
-        blank=gal_fits_bitpix_blank(bitpix);
+        blank=gal_fits_datatype_blank( gal_fits_bitpix_to_dtype(bitpix) );
         if(fits_write_key(fptr, datatype, "BLANK", blank,
                           "Pixels with no data.", &status) )
           gal_fits_io_error(status, "adding the BLANK keyword");
diff --git a/lib/gnuastro/fits.h b/lib/gnuastro/fits.h
index b061bf5..bb20701 100644
--- a/lib/gnuastro/fits.h
+++ b/lib/gnuastro/fits.h
@@ -26,18 +26,36 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <math.h>
 #include <float.h>
 #include <stdint.h>
+#include <limits.h>
 
 #include <fitsio.h>
 #include <wcslib/wcs.h>
 #include <wcslib/wcshdr.h>
 #include <wcslib/wcsfix.h>
+#include <gsl/gsl_complex.h>
+
+
+
+
+
+/* Order is based on the CFITSIO manual. Note that for the unsigned types
+   or small types (like char), the maximum value is considered as a blank
+   value, since the minimum value of an unsigned type is zero and zero is
+   often meaningful in contexts were unsigned values are used. */
+#define GAL_FITS_STRING_BLANK            NULL
+#define GAL_FITS_BYTE_BLANK              UCHAR_MAX
+#define GAL_FITS_LOGICAL_BLANK           SCHAR_MAX
+#define GAL_FITS_SHORT_BLANK             INT16_MIN
+#define GAL_FITS_LONG_BLANK              INT32_MIN
+#define GAL_FITS_LLONG_BLANK             INT64_MIN
+#define GAL_FITS_FLOAT_BLANK             NAN
+#define GAL_FITS_TCOMPLEX_BLANK(cptr)    GSL_SET_COMPLEX((cptr),NAN,NAN)
+#define GAL_FITS_TDBLCOMPLEX_BLANK(cptr) GSL_SET_COMPLEX((cptr),NAN,NAN)
+#define GAL_FITS_INT_BLANK               INT_MIN
+#define GAL_FITS_SBYTE_BLANK             SCHAR_MAX
+#define GAL_FITS_UINT_BLANK              UINT_MAX
+#define GAL_FITS_USHORT_BLANK            USHRT_MAX
 
-#define GAL_FITS_STRING_BLANK   NULL
-#define GAL_FITS_BYTE_BLANK     UCHAR_MAX /* 0 is often meaningful here! */
-#define GAL_FITS_SHORT_BLANK    INT16_MIN
-#define GAL_FITS_LONG_BLANK     INT32_MIN
-#define GAL_FITS_LLONG_BLANK    INT64_MIN
-#define GAL_FITS_FLOAT_BLANK    NAN
 
 
 
@@ -160,7 +178,7 @@ gal_fits_copyright_end(fitsfile *fptr,
  ******************        Read/Write        *****************
  *************************************************************/
 void *
-gal_fits_bitpix_blank(int bitpix);
+gal_fits_datatype_blank(int datatype);
 
 void
 gal_fits_convert_blank(void *array, int bitpix, size_t size, void *value);
diff --git a/src/imgcrop/crop.c b/src/imgcrop/crop.c
index 2d234fd..9ad5534 100644
--- a/src/imgcrop/crop.c
+++ b/src/imgcrop/crop.c
@@ -295,8 +295,8 @@ polygonmask(struct cropparams *crp, void *array, long 
*fpixel_i,
   LONGLONG *Lb, *La=array;
   unsigned char *bb, *ba=array;
   double *db, *ipolygon, point[2], *da=array;
-  int outpolygon=crp->p->outpolygon, bitpix=crp->p->bitpix;
   size_t i, *ordinds, size=s0*s1, nvertices=crp->p->nvertices;
+  int outpolygon=crp->p->outpolygon, datatype=crp->p->datatype;
 
 
   /* First of all, allocate enough space to put a copy of the input
@@ -325,10 +325,10 @@ polygonmask(struct cropparams *crp, void *array, long 
*fpixel_i,
 
   /* Go over all the pixels in the image and if they are within the
      polygon keep them if the user has asked for it.*/
-  switch(bitpix)
+  switch(datatype)
     {
-    case BYTE_IMG:
-      bb=gal_fits_bitpix_blank(bitpix);
+    case TBYTE:
+      bb=gal_fits_datatype_blank(datatype);
       for(i=0;i<size;++i)
         {
           point[0]=i%s1+1; point[1]=i/s1+1;
@@ -337,8 +337,8 @@ polygonmask(struct cropparams *crp, void *array, long 
*fpixel_i,
         }
       free(bb);
       break;
-    case SHORT_IMG:
-      sb=gal_fits_bitpix_blank(bitpix);
+    case TSHORT:
+      sb=gal_fits_datatype_blank(datatype);
       for(i=0;i<size;++i)
         {
           point[0]=i%s1+1; point[1]=i/s1+1;
@@ -347,8 +347,8 @@ polygonmask(struct cropparams *crp, void *array, long 
*fpixel_i,
         }
       free(sb);
       break;
-    case LONG_IMG:
-      lb=gal_fits_bitpix_blank(bitpix);
+    case TLONG:
+      lb=gal_fits_datatype_blank(datatype);
       for(i=0;i<size;++i)
         {
           point[0]=i%s1+1; point[1]=i/s1+1;
@@ -357,8 +357,8 @@ polygonmask(struct cropparams *crp, void *array, long 
*fpixel_i,
         }
       free(lb);
       break;
-    case LONGLONG_IMG:
-      Lb=gal_fits_bitpix_blank(bitpix);
+    case TLONGLONG:
+      Lb=gal_fits_datatype_blank(datatype);
       for(i=0;i<size;++i)
         {
           point[0]=i%s1+1; point[1]=i/s1+1;
@@ -367,8 +367,8 @@ polygonmask(struct cropparams *crp, void *array, long 
*fpixel_i,
         }
       free(Lb);
       break;
-    case FLOAT_IMG:
-      fb=gal_fits_bitpix_blank(bitpix);
+    case TFLOAT:
+      fb=gal_fits_datatype_blank(datatype);
       for(i=0;i<size;++i)
         {
           point[0]=i%s1+1; point[1]=i/s1+1;
@@ -377,8 +377,8 @@ polygonmask(struct cropparams *crp, void *array, long 
*fpixel_i,
         }
       free(fb);
       break;
-    case DOUBLE_IMG:
-      db=gal_fits_bitpix_blank(bitpix);
+    case TDOUBLE:
+      db=gal_fits_datatype_blank(datatype);
       for(i=0;i<size;++i)
         {
           point[0]=i%s1+1; point[1]=i/s1+1;
@@ -390,8 +390,8 @@ polygonmask(struct cropparams *crp, void *array, long 
*fpixel_i,
     default:
       error(EXIT_FAILURE, 0, "a bug! Please contact us at %s, so we "
             "can fix the problem. For some reason, an unrecognized "
-            "bitpix value (%d) has been seen in polygonmask (crop.c)",
-            PACKAGE_BUGREPORT, bitpix);
+            "datatype value (%d) has been seen in polygonmask (crop.c)",
+            PACKAGE_BUGREPORT, datatype);
     }
 
   /* Clean up: */
diff --git a/src/imgcrop/ui.c b/src/imgcrop/ui.c
index c10e2ac..6d2f642 100644
--- a/src/imgcrop/ui.c
+++ b/src/imgcrop/ui.c
@@ -648,7 +648,7 @@ preparearrays(struct imgcropparams *p)
         {
           firstbitpix=p->bitpix;
           p->datatype=gal_fits_bitpix_to_dtype(p->bitpix);
-          p->bitnul=gal_fits_bitpix_blank(p->bitpix);
+          p->bitnul=gal_fits_datatype_blank(p->datatype);
         }
       else if(firstbitpix!=p->bitpix)
         error(EXIT_FAILURE, 0, "%s: BITPIX=%d. Previous images had a "