[gnuastro-commits] master 41a52de 1/2: Library (wcs): New function to clean floating point errors in WCS

[Mohammad Akhlaghi]

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

    Library (wcs): New function to clean floating point errors in WCS
    
    Until now, extremely small values could be printed in the output WCS
    structure (for example PC elements on the scale of 1e-17, or extremely
    small differences in CDELT elements).
    
    With this commit, the new 'gal_wcs_clean_errors' function is added. By
    default, it is called just before writing the WCS by Gnuastro's FITS
    library. It will round such values to be equal in the CDELT and to be 0, 1
    or -1 in the PC matrix.
---
 NEWS               |  5 +++++
 doc/gnuastro.texi  | 11 +++++++++++
 lib/fits.c         |  3 +++
 lib/gnuastro/wcs.h |  5 +++++
 lib/wcs.c          | 46 +++++++++++++++++++++++++++++++++++++++++++++-
 5 files changed, 69 insertions(+), 1 deletion(-)

diff --git a/NEWS b/NEWS
index fb0a459..dc4cba3 100644
--- a/NEWS
+++ b/NEWS
@@ -14,6 +14,11 @@ See the end of the file for license conditions.
      will be read as floating point. This also applies to column arithmetic
      in Table.
 
+  Library:
+   - GAL_WCS_FLTERROR: Limit to identify a floating point error for WCS.
+   - gal_wcs_clean_errors: clean major WCS components from errors specified
+       by the FITS keyword 'CRDER', or floating point errors.
+
 ** Removed features
 
 ** Changed features
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index eb45d55..c3a034a 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -22105,6 +22105,10 @@ section are mainly just wrappers over CFITSIO, WCSLIB 
and GSL library
 functions to help in common applications.
 
 
+@deffn Macro GAL_WCS_FLTERROR
+Limit of rounding for floating point errors.
+@end deffn
+
 @deftypefun {struct wcsprm *} gal_wcs_read_fitsptr (fitsfile @code{*fptr}, 
size_t @code{hstartwcs}, size_t @code{hendwcs}, int @code{*nwcs})
 [@strong{Not thread-safe}] Return the WCSLIB @code{wcsprm} structure that
 is read from the CFITSIO @code{fptr} pointer to an opened FITS file. Also
@@ -22171,6 +22175,13 @@ several methods to store the matrix. The output is an 
allocated square
 matrix with each side equal to the number of dimensions.
 @end deftypefun
 
+@deftypefun void gal_wcs_clean_small_errors (struct wcsprm @code{*wcs})
+Errors can make small differences between the pixel-scale elements 
(@code{CDELT}) and can also lead to extremely small values in the @code{PC} 
matrix.
+With this function, such errors will be ``cleaned'' as follows: 1) if the 
maximum difference between the @code{CDELT} elements is smaller than the 
reference error, it will be set to the mean value.
+When the FITS keyword @code{CRDER} (optional) is defined it will be used as a 
reference, if not the default value is @code{GAL_WCS_FLTERROR}.
+2) If any of the PC elements differ from 0, 1 or -1 by less than 
@code{GAL_WCS_FLTERROR}, they will be rounded to the respective value.
+@end deftypefun
+
 @deftypefun void gal_wcs_decompose_pc_cdelt (struct wcsprm @code{*wcs})
 Decompose the @code{PCi_j} and @code{CDELTi} elements of
 @code{wcs}. According to the FITS standard, in the @code{PCi_j} WCS
diff --git a/lib/fits.c b/lib/fits.c
index 96682ba..bc71103 100644
--- a/lib/fits.c
+++ b/lib/fits.c
@@ -2180,6 +2180,9 @@ gal_fits_img_write_to_ptr(gal_data_t *input, char 
*filename)
       /* Decompose the 'PCi_j' matrix and 'CDELTi' vector. */
       gal_wcs_decompose_pc_cdelt(towrite->wcs);
 
+      /* Clean up small errors in the PC matrix and CDELT values. */
+      gal_wcs_clean_errors(towrite->wcs);
+
       /* Convert the WCS information to text. */
       status=wcshdo(WCSHDO_safe, towrite->wcs, &nkeyrec, &wcsstr);
       if(status)
diff --git a/lib/gnuastro/wcs.h b/lib/gnuastro/wcs.h
index 840e6d2..f06075e 100644
--- a/lib/gnuastro/wcs.h
+++ b/lib/gnuastro/wcs.h
@@ -30,6 +30,8 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 #include <gnuastro/data.h>
 
+/* Assumed floating point error in the WCS-related functionality. */
+#define GAL_WCS_FLTERROR 1e-12
 
 /* C++ Preparations */
 #undef __BEGIN_C_DECLS
@@ -83,6 +85,9 @@ double *
 gal_wcs_warp_matrix(struct wcsprm *wcs);
 
 void
+gal_wcs_clean_errors(struct wcsprm *wcs);
+
+void
 gal_wcs_decompose_pc_cdelt(struct wcsprm *wcs);
 
 double
diff --git a/lib/wcs.c b/lib/wcs.c
index 030ba83..703d2e8 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -32,6 +32,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <assert.h>
 
 #include <gsl/gsl_linalg.h>
+#include <wcslib/wcsmath.h>
 
 #include <gnuastro/wcs.h>
 #include <gnuastro/tile.h>
@@ -594,6 +595,46 @@ gal_wcs_warp_matrix(struct wcsprm *wcs)
 
 
 
+/* Clean up small/negligible errros that are clearly caused by measurement
+   errors in the PC and CDELT elements. */
+void
+gal_wcs_clean_errors(struct wcsprm *wcs)
+{
+  double crdcheck=NAN;
+  size_t i, crdnum=0, ndim=wcs->naxis;
+  double mean, crdsum=0, sum=0, min=FLT_MAX, max=0;
+  double *pc=wcs->pc, *cdelt=wcs->cdelt, *crder=wcs->crder;
+
+  /* First clean up CDELT: if the CRDER keyword is set, then we'll use that
+     as a reference, if not, we'll use the absolute floating point error
+     defined in 'GAL_WCS_FLTERR'. */
+  for(i=0; i<ndim; ++i)
+    {
+      sum+=cdelt[i];
+      if(cdelt[i]>max) max=cdelt[i];
+      if(cdelt[i]<min) min=cdelt[i];
+      if(crder[i]!=UNDEFINED) {++crdnum; crdsum=crder[i];}
+    }
+  mean=sum/ndim;
+  crdcheck = crdnum ? crdsum/crdnum : GAL_WCS_FLTERROR;
+  if( (max-min)/mean < crdcheck )
+    for(i=0; i<ndim; ++i)
+      cdelt[i]=mean;
+
+  /* Now clean up the PC elements. If the diagonal elements are too close
+     to 0, 1, or -1, set them to 0 or 1 or -1. */
+  if(pc)
+    for(i=0;i<ndim*ndim;++i)
+      {
+        if(      fabs(pc[i] -  0 ) < GAL_WCS_FLTERROR ) pc[i]=0;
+        else if( fabs(pc[i] -  1 ) < GAL_WCS_FLTERROR ) pc[i]=1;
+        else if( fabs(pc[i] - -1 ) < GAL_WCS_FLTERROR ) pc[i]=-1;
+      }
+}
+
+
+
+
 
 /* According to the FITS standard, in the 'PCi_j' WCS formalism, the matrix
    elements m_{ij} are encoded in the 'PCi_j' keywords and the scale
@@ -747,7 +788,10 @@ gal_wcs_pixel_scale(struct wcsprm *wcs)
           }
 
       /* Do the check, print warning and make correction. */
-      if(maxrow!=minrow && maxrow/minrow>1e5 && warning_printed==0)
+      if(maxrow!=minrow
+         && maxrow/minrow>1e5    /* The difference between elements is large */
+         && maxrow/minrow<GAL_WCS_FLTERROR
+         && warning_printed==0)
         {
           fprintf(stderr, "\nWARNING: The input WCS matrix (possibly taken "
                   "from the FITS header keywords starting with 'CD' or 'PC') "