[gnuastro-commits] master c9caf94: Library (fits.h): wrapper over WCSLIB bug (present from 2015)

[Mohammad Akhlaghi]

branch: master
commit c9caf944f12ca574f7ae095560a71b360604efbe
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>

    Library (fits.h): wrapper over WCSLIB bug (present from 2015)
    
    Until now, WCSLIB would wrongly write the CDELT values as '0.0' for some
    types of data-cubes (where the scales of the CDELT values in each dimension
    was large, but not large enough to warrant a different printing
    format). After contacting Mark Calabretta (author of WCSLIB), we discovered
    that this bug bug has been present in WCSLIB since its version 5.9
    (released on 2015/07/21). For more on the precise cause of the function see
    the comments above the new 'fits_bug_wrapper_cdelt_zero' function.
    
    Mark did mention that it will be fixed in the next release of
    WCSLIB. However, it will take a long time for many users (actually their
    package managers) to update. We therefore need to write a wrapper to fix
    the situation when it happens.
    
    With this commit, a bug-fix-wrapper function has been added just after
    writing the WCS keywords into the FITS file: upon writing any of the CDELT
    keyword lines into the FITS pointer, it will check if the written value,
    and the value inside the 'wcsprm' structure differ by more than 10
    decimals. When they do differ, by more than this, this function will use
    CFITSIO to update the keyword value in the key header.
    
    The operational footprint of this check is very small (only being activated
    on datasets that have more than 2 dimensions, and only doing the check on
    'CDELT' keywords, not the others).
    
    Also, since this is the first change in Gnuastro's library after the
    release of Gnuastro 0.14, I have incremented the Gnuastro library Soname
    version and added the top-level lines in the NEWS file for the next
    release.
---
 NEWS                         | 17 ++++++++
 configure.ac                 |  2 +-
 doc/announce-acknowledge.txt |  2 +
 doc/gnuastro.texi            | 11 +++--
 lib/fits.c                   | 95 ++++++++++++++++++++++++++++++++++++++++++--
 lib/gnuastro/fits.h          |  3 +-
 lib/wcs.c                    |  2 +-
 7 files changed, 120 insertions(+), 12 deletions(-)

diff --git a/NEWS b/NEWS
index ee96066..dd3b9a5 100644
--- a/NEWS
+++ b/NEWS
@@ -3,6 +3,23 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
 Copyright (C) 2015-2021 Free Software Foundation, Inc.
 See the end of the file for license conditions.
 
+* Noteworthy changes in release X.X (library X.X.X) (XXXX-XX-XX) [stable]
+
+** New features
+
+** Removed features
+
+** Changed features
+
+  Library:
+   - gal_fits_key_write_wcsstr: also takes WCS structure as argument.
+
+** Bugs fixed
+
+
+
+
+
 * Noteworthy changes in release 0.14 (library 12.0.0) (2021-01-25) [stable]
 
 ** New features
diff --git a/configure.ac b/configure.ac
index abeb585..1cc4d68 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=12
+GAL_CURRENT=13
 GAL_REVISION=0
 GAL_AGE=0
 GAL_LT_VERSION="${GAL_CURRENT}:${GAL_REVISION}:${GAL_AGE}"
diff --git a/doc/announce-acknowledge.txt b/doc/announce-acknowledge.txt
index 8622869..6aa0aca 100644
--- a/doc/announce-acknowledge.txt
+++ b/doc/announce-acknowledge.txt
@@ -1,5 +1,7 @@
 Alphabetically ordered list to acknowledge in the next release.
 
+Mark Calabretta
+
 
 
 
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 5064377..ee2dd1b 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -24027,12 +24027,11 @@ break it into several keywords (breaking up the 
string on directory
 separators).
 @end deftypefun
 
-@deftypefun void gal_fits_key_write_wcsstr (fitsfile @code{*fptr}, char 
@code{*wcsstr}, int @code{nkeyrec})
-Write the WCS header string (produced with WCSLIB's @code{wcshdo} function)
-into the CFITSIO @code{fitsfile} pointer. @code{nkeyrec} is the number of
-FITS header keywords in @code{wcsstr}. This function will put a few blank
-keyword lines along with a comment @code{WCS information} before writing
-each keyword record.
+@deftypefun void gal_fits_key_write_wcsstr (fitsfile @code{*fptr}, struct 
wcsprm @code{wcs}, char @code{*wcsstr}, int @code{nkeyrec})
+
+Write the WCS header string (produced with WCSLIB's @code{wcshdo} function) 
into the CFITSIO @code{fitsfile} pointer.
+@code{nkeyrec} is the number of FITS header keywords in @code{wcsstr}.
+This function will put a few blank keyword lines along with a comment 
@code{WCS information} before writing each keyword record.
 @end deftypefun
 
 @deftypefun void gal_fits_key_write (gal_fits_list_key_t @code{**keylist}, 
char @code{*title}, char @code{*filename}, char @code{*hdu})
diff --git a/lib/fits.c b/lib/fits.c
index a53a426..6dac672 100644
--- a/lib/fits.c
+++ b/lib/fits.c
@@ -1714,9 +1714,95 @@ gal_fits_key_write_filename(char *keynamebase, char 
*filename,
 
 
 
+/* A bug was found in WCSLIB that has existed since WCSLIB 5.9 (released on
+   2015/07/21) and will be fixed in the version after WCSLIB 7.3.1
+   (released in 2020). However, it will take time for many user package
+   managers to update their WCSLIB version. So we need to check if that bug
+   has occurred here and fix it manually.
+
+   In short the bug was originally seen on a dataset with this input CDELT
+   values:
+
+     CDELT1  =            0.0007778 / [deg]
+     CDELT2  =            0.0007778 / [deg]
+     CDELT3  =        30000.0000000 / [Hz]
+
+   The values are read into the 'wcsprm' structure properly, but upon
+   writing into the keyword string structure, the 'CDELT1' and 'CDELT2'
+   values are printed as 0. Mark Calabretta (creator of WCSLIB) described
+   the issue as follows:
+
+        """wcshdo() tries to find a single sprintf() floating point format
+        to use for groups of keywords, such as CDELTj as a group, or
+        CRPIXi, PCi_j, and CRVALj, each as separate groups.  It aims to
+        present the keyvalues in human-readable form, i.e. with decimal
+        points lined up, no unnecessary trailing zeroes, and avoiding
+        exponential ('E') format where its use is not warranted.
+
+        The problem here arose because of the large range of CDELT values
+        formatted using 'f' format, but not being so large that it would
+        force wcshdo() to revert to 'E' format.  There is also the
+        troubling possibility that in less extreme cases, precision of the
+        CDELT (or other) values could be lost without being noticed."""
+
+   To implement the check we will follow this logic: large dimensional
+   differences will not commonly happen in 2D datasets, so we will only
+   attempt the check in 3D datasets. We'll read each written CDELT value
+   with CFITSIO and if its zero, we'll correct it. */
+static void
+fits_bug_wrapper_cdelt_zero(fitsfile *fptr, struct wcsprm *wcs, char *keystr)
+{
+  size_t dim;
+  char *keyname;
+  double keyvalue;
+  int status=0, datatype=TDOUBLE;
+
+  /* Only do this check when we have more than two dimensions. */
+  if(wcs->naxis>2 && !strncmp(keystr, "CDELT", 5))
+    {
+      /* Find the dimension number and keyword string. This can later be
+         improved/generalized by actually parsing the keyword name to
+         extract the dimension, but I didn't have time to implement it in
+         the first implementation. It will rarely be necessary to go beyond
+         the third dimension, so this has almost no extra burden on the
+         computer's processing. */
+      if(      !strncmp(keystr, "CDELT1", 6) ) { keyname="CDELT1"; dim=0; }
+      else if( !strncmp(keystr, "CDELT2", 6) ) { keyname="CDELT2"; dim=1; }
+      else if( !strncmp(keystr, "CDELT3", 6) ) { keyname="CDELT3"; dim=2; }
+      else if( !strncmp(keystr, "CDELT4", 6) ) { keyname="CDELT4"; dim=3; }
+      else if( !strncmp(keystr, "CDELT5", 6) ) { keyname="CDELT5"; dim=4; }
+      else if( !strncmp(keystr, "CDELT6", 6) ) { keyname="CDELT6"; dim=5; }
+      else if( !strncmp(keystr, "CDELT7", 6) ) { keyname="CDELT7"; dim=6; }
+      else if( !strncmp(keystr, "CDELT8", 6) ) { keyname="CDELT8"; dim=7; }
+      else if( !strncmp(keystr, "CDELT9", 6) ) { keyname="CDELT9"; dim=8; }
+      else
+        error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix "
+              "the problem. There appears to be more than 9 dimensions in "
+              "the input dataset", __func__, PACKAGE_BUGREPORT);
+
+      /* Read the keyword value. */
+      fits_read_key(fptr, datatype, keyname, &keyvalue, NULL, &status);
+      gal_fits_io_error(status, NULL);
+
+      /* If the written value is not the same by more than 10 decimals,
+         re-write the value. */
+      if( fabs( wcs->cdelt[dim] - keyvalue ) > 1e-10  )
+        {
+          fits_update_key(fptr, datatype, keyname, &wcs->cdelt[dim],
+                          NULL, &status);
+          gal_fits_io_error(status, NULL);
+        }
+    }
+}
+
+
+
+
+
 /* Write the WCS header string into a FITS files*/
 void
-gal_fits_key_write_wcsstr(fitsfile *fptr, char *wcsstr, int nkeyrec)
+gal_fits_key_write_wcsstr(fitsfile *fptr, struct wcsprm *wcs,
+                          char *wcsstr, int nkeyrec)
 {
   size_t i;
   int status=0;
@@ -1738,7 +1824,10 @@ gal_fits_key_write_wcsstr(fitsfile *fptr, char *wcsstr, 
int nkeyrec)
          separate keyword along with all other important software, so it is
          redundant and just makes the keywrods hard to read by eye.*/
       if( keystart[0]!=' ' && strncmp(keystart, "COMMENT", 7) )
-        fits_write_record(fptr, keystart, &status);
+        {
+          fits_write_record(fptr, keystart, &status);
+          if(wcs) fits_bug_wrapper_cdelt_zero(fptr, wcs, keystart);
+        }
     }
   gal_fits_io_error(status, NULL);
 }
@@ -2500,7 +2589,7 @@ gal_fits_img_write_corr_wcs_str(gal_data_t *input, char 
*filename,
   fptr=gal_fits_img_write_to_ptr(input, filename);
 
   /* Write the WCS headers into the FITS file. */
-  gal_fits_key_write_wcsstr(fptr, wcsstr, nkeyrec);
+  gal_fits_key_write_wcsstr(fptr, NULL, wcsstr, nkeyrec);
 
   /* Update the CRPIX keywords. Note that we don't want to change the
      values in the WCS information of gal_data_t. Because, it often happens
diff --git a/lib/gnuastro/fits.h b/lib/gnuastro/fits.h
index 24c119b..aebd667 100644
--- a/lib/gnuastro/fits.h
+++ b/lib/gnuastro/fits.h
@@ -238,7 +238,8 @@ gal_fits_key_write_filename(char *keynamebase, char 
*filename,
                             gal_fits_list_key_t **list, int top1end0);
 
 void
-gal_fits_key_write_wcsstr(fitsfile *fptr, char *wcsstr, int nkeyrec);
+gal_fits_key_write_wcsstr(fitsfile *fptr, struct wcsprm *wcs,
+                          char *wcsstr, int nkeyrec);
 
 void
 gal_fits_key_write(gal_fits_list_key_t **keylist, char *title,
diff --git a/lib/wcs.c b/lib/wcs.c
index 028e2ab..0d2a508 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -517,7 +517,7 @@ gal_wcs_write_in_fitsptr(fitsfile *fptr, struct wcsprm *wcs)
           "wcshdu ERROR %d: %s", __func__, status, wcs_errmsg[status]);
   else
     {
-      gal_fits_key_write_wcsstr(fptr, wcsstr, nkeyrec);
+      gal_fits_key_write_wcsstr(fptr, wcs, wcsstr, nkeyrec);
       free(wcsstr);
     }
   status=0;