[gnuastro-commits] master 4b747b0: Check on floating point errors in pixel scale calculation

[Mohammad Akhlaghi]

branch: master
commit 4b747b037b280b1c5737bff227a3bb9a04375d21
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>

    Check on floating point errors in pixel scale calculation
    
    See the added error description for the full discussion: in short, some
    programs include floating point errors in their WCS FITS header keywords
    and these can bias the pixel scale measurements (and produce an error
    message sometimes). So with this commit, a warning will now be printed to
    guide the users on the (possible) problem and how to solve it.
    
    This was brought up in a discussion with Leindert Boogaard.
---
 THANKS                       |  1 +
 doc/announce-acknowledge.txt |  1 +
 lib/wcs.c                    | 59 ++++++++++++++++++++++++++++++++++++++++++--
 3 files changed, 59 insertions(+), 2 deletions(-)

diff --git a/THANKS b/THANKS
index 5f34189..1afb512 100644
--- a/THANKS
+++ b/THANKS
@@ -18,6 +18,7 @@ support in Gnuastro. The list is ordered alphabetically (by 
family name).
     Marjan Akbari                        address@hidden
     Karl Berry                           address@hidden
     Roland Bacon                         address@hidden
+    Leindert Boogaard                    address@hidden
     Nicolas Bouché                       address@hidden
     Fernando Buitrago                    address@hidden
     Adrian Bunk                          address@hidden
diff --git a/doc/announce-acknowledge.txt b/doc/announce-acknowledge.txt
index 880c69c..4e2bfea 100644
--- a/doc/announce-acknowledge.txt
+++ b/doc/announce-acknowledge.txt
@@ -1,5 +1,6 @@
 This file is meant to keep the names of the people who's help must be
 acknowledged in the next release.
 
+Leindert Boogaard
 Takashi Ichikawa
 David Valls-Gabaud
diff --git a/lib/wcs.c b/lib/wcs.c
index 1fe18da..5b438eb 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -438,10 +438,10 @@ double *
 gal_wcs_pixel_scale(struct wcsprm *wcs)
 {
   gsl_vector S;
-  double *a, *out;
-  int permute_set;
   gsl_matrix A, V;
   size_t i, j, n=wcs->naxis;
+  double *a, *out, maxrow, minrow;
+  int permute_set, warning_printed;
   double *v=gal_data_malloc_array(GAL_TYPE_FLOAT64, n*n, __func__, "v");
   size_t *permutation=gal_data_malloc_array(GAL_TYPE_SIZE_T, n, __func__,
                                             "permutation");
@@ -455,6 +455,61 @@ gal_wcs_pixel_scale(struct wcsprm *wcs)
   a=gal_wcs_warp_matrix(wcs);
 
 
+  /* To avoid confusing issues with floating point errors being written in
+     the non-diagonal elements of the FITS header PC or CD matrices, we
+     need to check if the minimum and maximum values in each row are not
+     several orders of magnitude apart.
+
+     Note that in some cases (for example a spectrum), one axis might be in
+     degrees (value around 1e-5) and the other in angestroms (value around
+     1e-10). So we can't look at the minimum and maximum of the whole
+     matrix. However, in such cases, people will probably not warp/rotate
+     the image to mix the coordinates. So the important thing to check is
+     the minimum and maximum (non-zero) values in each row. */
+  warning_printed=0;
+  for(i=0;i<n;++i)
+    {
+      /* Find the minimum and maximum values in each row. */
+      minrow=FLT_MAX;
+      maxrow=-FLT_MAX;
+      for(j=0;j<n;++j)
+        if(a[i*n+j]!=0.0) /* We aren't concerned with 0 valued elements. */
+          {
+            /* We have to use absolutes because in cases like RA, the
+               diagonal values in different rows may have different signs*/
+            if(fabs(a[i*n+j])<minrow) minrow=fabs(a[i*n+j]);
+            if(fabs(a[i*n+j])>maxrow) maxrow=fabs(a[i*n+j]);
+          }
+
+      /* Do the check, print warning and make correction. */
+      if(maxrow!=minrow && maxrow/minrow>1e4 && warning_printed==0)
+        {
+          fprintf(stderr, "\nWARNING: The input WCS matrix (possibly taken "
+                  "from the FITS header keywords starting with `CD' or `PC') "
+                  "contains values with very different scales (more than "
+                  "10^4 different). This is probably due to floating point "
+                  "errors. These values might bias the pixel scale (and "
+                  "subsequent) calculations.\n\n"
+                  "You can see the respective matrix with one of the "
+                  "following two commands (depending on how the FITS file "
+                  "was written). Recall that if the desired extension/HDU "
+                  "isn't the default, you can choose it with the `--hdu' "
+                  "(or `-h') option before the `|' sign in these commands."
+                  "\n\n"
+                  "    $ astfits file.fits -p | grep 'PC._.'\n"
+                  "    $ astfits file.fits -p | grep 'CD._.'\n\n"
+                  "You can delete the ones with obvious floating point "
+                  "error values using the following command (assuming you "
+                  "want to delete `CD1_2' and `CD2_1') and rerun the your "
+                  "command to remove this warning message and possibly "
+                  "correct errors that might be reported due to it.\n\n"
+                  "    $ astfits file.fits --delete=CD1_2 --delete=CD2_1\n\n"
+                  );
+          warning_printed=1;
+        }
+    }
+
+
   /* Fill in the necessary GSL vector and Matrix structures. */
   S.size=n;     S.stride=1;                S.data=pixscale->array;
   V.size1=n;    V.size2=n;    V.tda=n;     V.data=v;