[gnuastro-commits] master 36f9ad2: Integer blank values in arithmetic and the FITS BLANK keyword

[Mohammad Akhlaghi]

branch: master
commit 36f9ad2d36192f55cd1e1319da8d4b38b2f47c33
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>

    Integer blank values in arithmetic and the FITS BLANK keyword
    
    This commit solves some issues related to integer blank values in the
    Arithmetic library and the BLANK keyword that must be written to the FITS
    headers. They were both found while trying to do a simple arithmetic
    operation on an unsigned 16-bit image and noticing that blank values are
    not being considered properly in both stages.
    
    For reading and writing a FITS image, we now have the
    `gal_fits_key_img_blank' function which will allocate and fill the memory
    for Gnuastro's blank values of any type as required by the FITS standard
    (because the BLANK value must be the value before BZERO and BSCALE are
    applied).
    
    In the process, I also noticed the integer limits reported in the manual
    for 8-bit and 16-bit signed integers needed corrections which are now
    applied.
    
    Also, in preparation for version 0.5, the library version has been
    incremented to 3.0.0.
    
    This fixes bug #52014.
---
 NEWS                    |  8 ++++++-
 bin/crop/ui.c           |  2 +-
 configure.ac            |  2 +-
 doc/gnuastro.texi       | 18 +++++++++++++--
 lib/arithmetic-binary.c | 50 ++++++++++++++++++++++++++++++++++-------
 lib/fits.c              | 59 ++++++++++++++++++++++++++++++++++++++++++++++++-
 lib/gnuastro/fits.h     |  3 +++
 7 files changed, 128 insertions(+), 14 deletions(-)

diff --git a/NEWS b/NEWS
index f58fc93..e3b992d 100644
--- a/NEWS
+++ b/NEWS
@@ -1,10 +1,13 @@
 GNU Astronomy Utilities NEWS                          -*- outline -*-
 
 
-* Noteworthy changes in release 0.4.XXX (library 2.0.0) (YYYY-MM-DD) [stable]
+* Noteworthy changes in release 0.4.XXX (library 3.0.0) (YYYY-MM-DD) [stable]
 
 ** New features
 
+  `gal_fits_key_img_blank': returns the value that must be used in the
+  BLANK keyword for the given type as defined by the FITS standard.
+
 ** Removed features
 
 ** Changed features
@@ -13,6 +16,9 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
 
   ConvertType crash when changing values (bug #52010).
 
+  Arithmetic not accounting for integer blank pixels in binary operators
+  (bug #52014).
+
 
 
 
diff --git a/bin/crop/ui.c b/bin/crop/ui.c
index 089a0d9..9c72b7a 100644
--- a/bin/crop/ui.c
+++ b/bin/crop/ui.c
@@ -882,7 +882,7 @@ ui_preparations(struct cropparams *p)
           /* Set the basic information. */
           firsttype = p->type;
           firstndim = img->ndim;
-          p->bitnul = gal_blank_alloc_write(p->type);
+          p->bitnul = gal_fits_key_img_blank(p->type);
 
           /* Make sure the number of dimensions is supported. */
           if(firstndim>MAXDIM)
diff --git a/configure.ac b/configure.ac
index 4b9d39f..73d8679 100644
--- a/configure.ac
+++ b/configure.ac
@@ -50,7 +50,7 @@ AC_CONFIG_MACRO_DIRS([bootstrapped/m4])
 
 # Library version, see the GNU Libtool manual ("Library interface versions"
 # section for the exact definition of each) for
-GAL_CURRENT=2
+GAL_CURRENT=3
 GAL_REVISION=0
 GAL_AGE=0
 GAL_LT_VERSION="${GAL_CURRENT}:${GAL_REVISION}:${GAL_AGE}"
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index f43ef19..cdde8ae 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -5437,7 +5437,7 @@ are inclusive.
 @item i8
 @itemx int8
 8-bit signed integers, range:@*
address@hidden to\ }2^7-1]} or @mymath{[-127\rm{\ to\ }127]}.
address@hidden to\ }2^7-1]} or @mymath{[-128\rm{\ to\ }127]}.
 
 @item u16
 @itemx uint16
@@ -5447,7 +5447,7 @@ are inclusive.
 @item i16
 @itemx int16
 16-bit signed integers, range:@* @mymath{[-2^{15}\rm{\ to\ }2^{15}-1]} or
address@hidden to\ }32768]}.
address@hidden to\ }32767]}.
 
 @item u32
 @itemx uint32
@@ -19409,6 +19409,20 @@ typedef struct gal_fits_list_key_t
 @end example
 @end deftp
 
address@hidden {void *} gal_fits_key_img_blank (uint8_t @code{type})
+Returns a pointer to an allocated space for the FITS @code{BLANK} header
+keyword when the input array has a type of @code{type}.
+
+According to the FITS standard: ``If the @code{BSCALE} and @code{BZERO}
+keywords do not have the default values of 1.0 and 0.0, respectively, then
+the value of the @code{BLANK} keyword must equal the actual value in the
+FITS data array that is used to represent an undefined pixel and not the
+corresponding physical value''. Therefore a special @code{BLANK} value is
+needed for datasets containing signed 8-bit integers and unsigned 16-bit,
+32-bit and 64-bit integers (types that are defined with @code{BSCALE} and
address@hidden in the FITS standard).
address@hidden deftypefun
+
 @deftypefun void gal_fits_key_clean_str_value (char @code{*string})
 Remove the single quotes and possible extra spaces around the keyword values
 that CFITSIO returns when reading a string keyword. CFITSIO doesn't remove
diff --git a/lib/arithmetic-binary.c b/lib/arithmetic-binary.c
index 02d7ce1..c5facc8 100644
--- a/lib/arithmetic-binary.c
+++ b/lib/arithmetic-binary.c
@@ -27,6 +27,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <error.h>
 #include <stdlib.h>
 
+#include <gnuastro/blank.h>
 #include <gnuastro/arithmetic.h>
 
 #include <gnuastro-internal/arithmetic-internal.h>
@@ -196,20 +197,41 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 /*************              High level macros           *****************/
 /************************************************************************/
 /* Final step to be used by all operators and all types. */
-#define BINARY_OP_OT_RT_LT_SET(OP, OT, RT, LT) {                   \
-    LT *la=l->array;                                               \
-    RT *ra=r->array;                                               \
-    OT *oa=o->array, *of=oa + o->size;                             \
-    if(l->size==r->size) do *oa = *la++ OP *ra++; while(++oa<of);  \
-    else if(l->size==1)  do *oa = *la   OP *ra++; while(++oa<of);  \
-    else                 do *oa = *la++ OP *ra;   while(++oa<of);  \
+#define BINARY_OP_OT_RT_LT_SET(OP, OT, RT, LT) {                        \
+    LT lb, *la=l->array;                                                \
+    RT rb, *ra=r->array;                                                \
+    OT ob, *oa=o->array, *of=oa + o->size;                              \
+    if(checkblank)                                                      \
+      {                                                                 \
+        gal_blank_write(&lb, l->type);                                  \
+        gal_blank_write(&rb, r->type);                                  \
+        gal_blank_write(&ob, o->type);                                  \
+        do                                                              \
+          {                                                             \
+            if(lb==lb && rb==rb)/* Both are integers.                */ \
+              *oa = (*la!=lb  && *ra!=rb)  ? *la OP *ra : ob ;          \
+            else if(lb==lb)     /* Only left operand is an integer.  */ \
+              *oa = (*la!=lb  && *ra==*ra) ? *la OP *ra : ob;           \
+            else                /* Only right operand is an integer. */ \
+              *oa = (*la==*la && *ra!=rb)  ? *la OP *ra : ob;           \
+            if(l->size>1) ++la;                                         \
+            if(r->size>1) ++ra;                                         \
+          }                                                             \
+        while(++oa<of);                                                 \
+      }                                                                 \
+    else                                                                \
+      {                                                                 \
+        if(l->size==r->size) do *oa = *la++ OP *ra++; while(++oa<of);   \
+        else if(l->size==1)  do *oa = *la   OP *ra++; while(++oa<of);   \
+        else                 do *oa = *la++ OP *ra;   while(++oa<of);   \
+      }                                                                 \
   }
 
 
 
 
 /* This is for operators like `&&' and `||', where the right operator is
-   not necessarily read (and thus incremented) incremented. */
+   not necessarily read (and thus incremented). */
 #define BINARY_OP_INCR_OT_RT_LT_SET(OP, OT, RT, LT) {                   \
     LT *la=l->array;                                                    \
     RT *ra=r->array;                                                    \
@@ -332,6 +354,7 @@ arithmetic_binary(int operator, uint8_t flags, gal_data_t 
*lo,
   /* Read the variable arguments. `lo' and `ro' keep the original data, in
      case their type isn't built (based on configure options are configure
      time). */
+  int checkblank;
   int32_t otype, final_otype;
   size_t out_size, minmapsize;
   gal_data_t *l, *r, *o=NULL, *tmp_o;
@@ -398,6 +421,17 @@ arithmetic_binary(int operator, uint8_t flags, gal_data_t 
*lo,
                        0, minmapsize, NULL, NULL, NULL );
 
 
+  /* See if we should check for blanks. When both types are floats, blanks
+     don't need to be checked (the floating point standard will do the job
+     for us). It is also not necessary to check blanks in bitwise
+     operators, but bitwise operators have their own macro
+     (`BINARY_OP_INCR_OT_RT_LT_SET') which doesn' use `checkblanks'.*/
+  checkblank = ((((l->type!=GAL_TYPE_FLOAT32    && l->type!=GAL_TYPE_FLOAT64)
+                  || (r->type!=GAL_TYPE_FLOAT32 && r->type!=GAL_TYPE_FLOAT64))
+                 && (gal_blank_present(l, 1) || gal_blank_present(r, 1)))
+                ? 1 : 0 );
+
+
   /* Start setting the operator and operands. */
   switch(l->type)
     {
diff --git a/lib/fits.c b/lib/fits.c
index 9df4228..b0f901e 100644
--- a/lib/fits.c
+++ b/lib/fits.c
@@ -794,6 +794,63 @@ gal_fits_hdu_open_format(char *filename, char *hdu, int 
img0_tab1)
 /**************************************************************/
 /**********            Header keywords             ************/
 /**************************************************************/
+/* Return allocated pointer to the blank value to use in a FITS file header
+   keyword. */
+void *
+gal_fits_key_img_blank(uint8_t type)
+{
+  void *out=NULL, *tocopy=NULL;
+  uint8_t u8=0;          /* Equivalent of minimum in signed   with BZERO. */
+  int16_t s16=INT16_MAX; /* Equivalend of maximum in unsigned with BZERO. */
+  int32_t s32=INT32_MAX; /* Equivalend of maximum in unsigned with BZERO. */
+  int64_t s64=INT64_MAX; /* Equivalend of maximum in unsigned with BZERO. */
+
+  switch(type)
+    {
+    /* Types with no special treatment: */
+    case GAL_TYPE_BIT:
+    case GAL_TYPE_UINT8:
+    case GAL_TYPE_INT16:
+    case GAL_TYPE_INT32:
+    case GAL_TYPE_INT64:
+    case GAL_TYPE_FLOAT32:
+    case GAL_TYPE_FLOAT64:
+    case GAL_TYPE_COMPLEX32:
+    case GAL_TYPE_COMPLEX64:
+    case GAL_TYPE_STRING:
+    case GAL_TYPE_STRLL:
+      out = gal_blank_alloc_write(type);
+      break;
+
+    /* Types that need values from their opposite-signed types. */
+    case GAL_TYPE_INT8:      tocopy=&u8;      break;
+    case GAL_TYPE_UINT16:    tocopy=&s16;     break;
+    case GAL_TYPE_UINT32:    tocopy=&s32;     break;
+    case GAL_TYPE_UINT64:    tocopy=&s64;     break;
+
+    default:
+      error(EXIT_FAILURE, 0, "%s: type code %u not recognized as a Gnuastro "
+            "data type", __func__, type);
+    }
+
+  /* If `gal_blank_alloc_write' wasn't used (copy!=NULL), then allocate the
+     necessary space and fill it in. Note that the width of the signed and
+     unsigned values doesn't differ, so we can use the actual input
+     type. */
+  if(tocopy)
+    {
+      out = gal_data_malloc_array(type, 1, __func__, "out");
+      memcpy(out, tocopy, gal_type_sizeof(type));
+    }
+
+  /* Return. */
+  return out;
+}
+
+
+
+
+
 /* CFITSIO doesn't remove the two single quotes around the string value, so
    the strings it reads are like: 'value ', or 'some_very_long_value'. To
    use the value, it is commonly necessary to remove the single quotes (and
@@ -1710,7 +1767,7 @@ gal_fits_img_write_to_ptr(gal_data_t *input, char 
*filename)
         break;
 
       default:
-        blank=gal_blank_alloc_write(towrite->type);
+        blank=gal_fits_key_img_blank(towrite->type);
         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 a44a3d5..2744062 100644
--- a/lib/gnuastro/fits.h
+++ b/lib/gnuastro/fits.h
@@ -159,6 +159,9 @@ gal_fits_hdu_open_format(char *filename, char *hdu, int 
img0_tab1);
 /**************************************************************/
 /**********            Header keywords             ************/
 /**************************************************************/
+void *
+gal_fits_key_img_blank(uint8_t type);
+
 void
 gal_fits_key_clean_str_value(char *string);