[gnuastro-commits] master 9c7d9a0 10/16: Blank replacement function now accepts datatype

[Mohammad Akhlaghi]

branch: master
commit 9c7d9a02d98a03ac2a703cd537683e013a0e18de
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>

    Blank replacement function now accepts datatype
    
    Following the discussion in commit 2071999 (Blank CFITSIO pointers set by
    datatype, not bitpix), the old `blanktovalue' function now also accepts all
    the CFITSIO dataypes. It is thus now named `gal_fits_blank_to_value' and
    included in the header file incase some programs might need it. It is
    currently not used anywhere, but can be useful so the effort that went to
    into its initial definition is now kept for future use, especially after
    libraries are installable.
    
    While working on this, the following two issues were also addressed in this
    commit:
    
      - The unsigned long datatype was missed in the commits since commit
        2071999, so it is now included in the functions dealing with type.
    
      - The CFITSIO specific Bitpix values are now also converted to datatype
        in `gal_fits_bitpix_to_dtype'.
    
      - Until now a single `GAL_FITS_FLOAT_BLANK' was used (assuming it is a
        value of NaN. But for some reason, it might be useful to have a
        different blank value for floats or doubles. Therefore a new
        `GAL_FITS_DOUBLE_BLANK' was also defined to be more general.
---
 lib/fits.c          |  196 +++++++++++++++++++++++++++++++++++++++++----------
 lib/gnuastro/fits.h |   29 ++++----
 2 files changed, 174 insertions(+), 51 deletions(-)

diff --git a/lib/fits.c b/lib/fits.c
index edee58d..e1343e6 100644
--- a/lib/fits.c
+++ b/lib/fits.c
@@ -167,6 +167,12 @@ gal_fits_bitpix_to_dtype(int bitpix)
       return TFLOAT;
     case DOUBLE_IMG:
       return TDOUBLE;
+    case SBYTE_IMG:
+      return TSBYTE;
+    case USHORT_IMG:
+      return TUSHORT;
+    case ULONG_IMG:
+      return TULONG;
     default:
       error(EXIT_FAILURE, 0, "bitpix value of %d not recognized",
             bitpix);
@@ -243,11 +249,12 @@ gal_fits_datatype_blank(int datatype)
   LONGLONG *L;
   float *f;
   double *d;
-  gsl_complex_float *complex;
-  gsl_complex *dblcomplex;
+  gsl_complex_float *cx;
+  gsl_complex *dcx;
   int *i;
   unsigned int *ui;
   unsigned short *us;
+  unsigned long *ul;
 
   errno=0;
   switch(datatype)
@@ -320,24 +327,24 @@ gal_fits_datatype_blank(int datatype)
       if(d==NULL)
         error(EXIT_FAILURE, errno, "%lu bytes for blank TDOUBLE",
               sizeof *d);
-      *d=GAL_FITS_FLOAT_BLANK;
+      *d=GAL_FITS_DOUBLE_BLANK;
       return d;
 
     case TCOMPLEX:
-      complex=malloc(sizeof *complex);
-      if(complex==NULL)
+      cx=malloc(sizeof *cx);
+      if(cx==NULL)
         error(EXIT_FAILURE, errno, "%lu bytes for blank TCOMPLEX",
-              sizeof *complex);
-      GAL_FITS_TCOMPLEX_BLANK(complex);
-      return complex;
+              sizeof *cx);
+      GSL_SET_COMPLEX(cx,GAL_FITS_FLOAT_BLANK,GAL_FITS_FLOAT_BLANK);
+      return cx;
 
     case TDBLCOMPLEX:
-      dblcomplex=malloc(sizeof *dblcomplex);
-      if(dblcomplex==NULL)
+      dcx=malloc(sizeof *dcx);
+      if(dcx==NULL)
         error(EXIT_FAILURE, errno, "%lu bytes for blank TDBLCOMPLEX",
-              sizeof *dblcomplex);
-      GAL_FITS_TCOMPLEX_BLANK(dblcomplex);
-      return dblcomplex;
+              sizeof *dcx);
+      GSL_SET_COMPLEX(dcx,GAL_FITS_DOUBLE_BLANK,GAL_FITS_DOUBLE_BLANK);
+      return dcx;
 
     case TINT:
       i=malloc(sizeof *i);
@@ -363,6 +370,14 @@ gal_fits_datatype_blank(int datatype)
       *ui=GAL_FITS_USHORT_BLANK;
       return us;
 
+    case TULONG:
+      ul=malloc(sizeof *ul);
+      if(ul==NULL)
+        error(EXIT_FAILURE, errno, "%lu bytes for blank TULONG",
+              sizeof *ul);
+      *ul=GAL_FITS_ULONG_BLANK;
+      return ul;
+
 
     default:
       error(EXIT_FAILURE, 0, "datatype value of %d not recognized",
@@ -464,6 +479,10 @@ gal_fits_datatype_alloc(size_t size, int datatype)
       size *= sizeof (unsigned short);
       break;
 
+    case TULONG:
+      size *= sizeof (unsigned long);
+      break;
+
     default:
       error(EXIT_FAILURE, 0, "datatype value of %d not recognized in "
             "gal_fits_datatype_alloc", datatype);
@@ -483,54 +502,155 @@ gal_fits_datatype_alloc(size_t size, int datatype)
 
 
 void
-blanktovalue(void *array, int bitpix, size_t size, void *value)
+gal_fits_blank_to_value(void *array, int datatype, size_t size, void *value)
 {
-  /* 'value' will only be read from one of these based on bitpix. Which the
-     caller assigned. If there is any problem, it is their responsability,
-     not this functions :-).*/
-  unsigned char *b, *bf, bv=*(uint8_t *) value;
-  short *s, *sf, sv=*(int16_t *) value;
-  long *l, *lf, lv=*(int32_t *) value;
-  LONGLONG *L, *Lf, Lv=*(int64_t *) value;
-  float   *f, *ff, fv=*(float   *) value;
-  double  *d, *df, dv=*(double  *) value;
+  /* 'value' will only be read from one of these based on the
+     datatype. Which the caller assigned. If there is any problem, it is
+     their responsability, not this function's.*/
+  unsigned char *b, *bf,        bv = *(uint8_t *) value;
+  char *c, *cf,                 cv = *(char *) value;
+  char **str, **strf,        *strv = *(char **) value;
+  short *s, *sf,                sv = *(int16_t *) value;
+  long *l, *lf,                 lv = *(int32_t *) value;
+  LONGLONG *L, *Lf,             Lv = *(int64_t *) value;
+  float   *f, *ff,              fv = *(float *) value;
+  double  *d, *df,              dv = *(double *) value;
+  gsl_complex_float *cx, *cxf, cxv = *(gsl_complex_float *) value;
+  gsl_complex *dcx, *dcxf,    dcxv = *(gsl_complex *) value;
+  int *in, *inf,               inv = *(int *) value;
+  unsigned int *ui, *uif,      uiv = *(unsigned int *) value;
+  unsigned short *us, *usf,    usv = *(unsigned short *) value;
+  unsigned long *ul, *ulf,     ulv = *(unsigned long *) value;
 
-  switch(bitpix)
+  switch(datatype)
     {
-    case BYTE_IMG:
+    case TBIT:
+      error(EXIT_FAILURE, 0, "Currently Gnuastro doesn't support TBIT "
+            "datatype, please get in touch with us to implement it.");
+
+    case TBYTE:
       bf=(b=array)+size;
-      do if(*b==GAL_FITS_BYTE_BLANK) *b=bv; while(++b<bf);
+      do if(*b==GAL_FITS_BYTE_BLANK) *b++=bv; while(b<bf);
       break;
 
-    case SHORT_IMG:
+
+    case TLOGICAL: case TSBYTE:
+      cf=(c=array)+size;
+      do if(*c==GAL_FITS_LOGICAL_BLANK) *c++=cv; while(c<cf);
+      break;
+
+
+    case TSTRING:
+      strf=(str=array)+size;
+      do if(*str==GAL_FITS_STRING_BLANK) *str++=strv; while(str<strf);
+      break;
+
+
+    case TSHORT:
       sf=(s=array)+size;
-      do if(*s==GAL_FITS_SHORT_BLANK) *s=sv; while(++s<sf);
+      do if(*s==GAL_FITS_SHORT_BLANK) *s++=sv; while(s<sf);
       break;
 
-    case LONG_IMG:
+
+    case TLONG:
       lf=(l=array)+size;
-      do if(*l==GAL_FITS_LONG_BLANK) *l=lv; while(++l<lf);
+      do if(*l==GAL_FITS_LONG_BLANK) *l++=lv; while(l<lf);
       break;
 
-    case LONGLONG_IMG:
+
+    case TLONGLONG:
       Lf=(L=array)+size;
-      do if(*L==GAL_FITS_LLONG_BLANK) *L=Lv; while(++L<Lf);
+      do if(*L==GAL_FITS_LLONG_BLANK) *L++=Lv; while(L<Lf);
       break;
 
-    case FLOAT_IMG:
+
+      /* Note that a NaN value is not equal to another NaN value, so we
+         can't use the easy check for cases were the blank value is
+         NaN. Also note that `isnan' is actually a macro, so it works for
+         both float and double types.*/
+    case TFLOAT:
       ff=(f=array)+size;
-      do if(*f==GAL_FITS_FLOAT_BLANK) *f=fv; while(++f<ff);
+      if(isnan(GAL_FITS_FLOAT_BLANK))
+        do if(isnan(*f)) *f++=fv; while(f<ff);
+      else
+        do if(*f==GAL_FITS_FLOAT_BLANK) *f++=fv; while(f<ff);
       break;
 
-    case DOUBLE_IMG:
+
+    case TDOUBLE:
       df=(d=array)+size;
-      do if(*d==GAL_FITS_FLOAT_BLANK) *d=dv; while(++d<df);
+      if(isnan(GAL_FITS_DOUBLE_BLANK))
+        do if(isnan(*d)) *d++=dv; while(d<df);
+      else
+        do if(*d==GAL_FITS_FLOAT_BLANK) *d++=dv; while(d<df);
+      break;
+
+
+    case TCOMPLEX:
+      cxf=(cx=array)+size;
+      if(isnan(GAL_FITS_FLOAT_BLANK))
+          do
+            if(isnan(GSL_COMPLEX_P_REAL(cx))
+               && isnan(GSL_COMPLEX_P_IMAG(cx)) )
+              GSL_SET_COMPLEX(cx, GSL_COMPLEX_P_REAL(&cxv),
+                              GSL_COMPLEX_P_IMAG(&cxv));
+          while(++cx<cxf);
+      else
+        do
+          if( GSL_COMPLEX_P_REAL(cx) == GAL_FITS_FLOAT_BLANK
+              && GSL_COMPLEX_P_IMAG(cx) == GAL_FITS_FLOAT_BLANK)
+            GSL_SET_COMPLEX(cx, GSL_COMPLEX_P_REAL(&cxv),
+                            GSL_COMPLEX_P_IMAG(&cxv));
+        while(++cx<cxf);
+      break;
+
+
+    case TDBLCOMPLEX:
+      dcxf=(dcx=array)+size;
+      if(isnan(GAL_FITS_DOUBLE_BLANK))
+          do
+            if(isnan(GSL_COMPLEX_P_REAL(dcx))
+               && isnan(GSL_COMPLEX_P_IMAG(dcx)) )
+              GSL_SET_COMPLEX(dcx, GSL_COMPLEX_P_REAL(&dcxv),
+                              GSL_COMPLEX_P_IMAG(&dcxv));
+          while(++dcx<dcxf);
+      else
+        do
+          if( GSL_COMPLEX_P_REAL(dcx) == GAL_FITS_FLOAT_BLANK
+              && GSL_COMPLEX_P_IMAG(dcx) == GAL_FITS_FLOAT_BLANK)
+            GSL_SET_COMPLEX(dcx, GSL_COMPLEX_P_REAL(&dcxv),
+                            GSL_COMPLEX_P_IMAG(&dcxv));
+        while(++dcx<dcxf);
+      break;
+
+
+    case TINT:
+      inf=(in=array)+size;
+      do if(*in==GAL_FITS_INT_BLANK) *in++=inv; while(in<inf);
+      break;
+
+
+    case TUINT:
+      uif=(ui=array)+size;
+      do if(*ui==GAL_FITS_UINT_BLANK) *ui++=uiv; while(ui<uif);
+      break;
+
+
+    case TUSHORT:
+      usf=(us=array)+size;
+      do if(*us==GAL_FITS_USHORT_BLANK) *us++=usv; while(us<usf);
+      break;
+
+
+    case TULONG:
+      ulf=(ul=array)+size;
+      do if(*ul==GAL_FITS_ULONG_BLANK) *ul++=ulv; while(ul<ulf);
       break;
 
     default:
-      error(EXIT_FAILURE, 0, "a bug! Bitpix value of %d not recognized. "
+      error(EXIT_FAILURE, 0, "a bug! datatype value of %d not recognized. "
             "This should not happen here (blanktovalue in fitsarrayvv.c). "
-            "Please contact us at %s to see how this happened", bitpix,
+            "Please contact us at %s to see how this happened", datatype,
             PACKAGE_BUGREPORT);
     }
 }
diff --git a/lib/gnuastro/fits.h b/lib/gnuastro/fits.h
index 8822b86..9e759c1 100644
--- a/lib/gnuastro/fits.h
+++ b/lib/gnuastro/fits.h
@@ -42,19 +42,19 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
    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
+#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_DOUBLE_BLANK     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_ULONG_BLANK      ULONG_MAX
 
 
 
@@ -183,6 +183,9 @@ gal_fits_datatype_blank(int datatype);
 void
 gal_fits_convert_blank(void *array, int bitpix, size_t size, void *value);
 
+void
+gal_fits_blank_to_value(void *array, int datatype, size_t size, void *value);
+
 int
 gal_fits_bitpix_to_dtype(int bitpix);