[gnuastro-commits] master d3deff2: Library (wcs.h): writing CD matrix when TPV distortion present

[Mohammad Akhlaghi]

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

    Library (wcs.h): writing CD matrix when TPV distortion present
    
    Until now, Gnuastro would simply follow WCSLIB's default output of writing
    the PC+CDELT matrices into the header keywords. However, we recently found
    out that the PC+CDELT matrix (when CDELT is not equal to 1) will cause
    problems in WCSLIB's TPV conversion. After contacting Mark Calabretta, we
    found out that this is a missing issue in WCSLIB and that it will be fixed
    in the next version of WCSLIB (7.4, to be released).
    
    This is the full reply of Mark: "the problem was that TPV assumes that
    CDi_ja is used to define the linear transformation. (Consequently, the
    independent variables of the TPV polynomial are intermediate world
    coordinates rather than intermediate pixel coordinates. That's the wrong
    way to do it, but can be excused as an historical anomaly.) Thus, unlike
    other sequent distortions, if PCi_ja is used instead of CDi_ja, then
    CDELTja must be applied before the TPV polynomial. WCSLIB wasn't doing
    that".
    
    Even though it will hopefully be fixed in WCSLIB 7.4, it will take a long
    time for the community to transition to it (maybe +3 years or so!). So it
    is safer if we also implement a correction in Gnuastro and avoid problems
    for users who have older WCSLIB versions.
    
    With this commit, if a TPV distortion is present, we'll multiply the CDELT
    matrix with the PC matrix (so effectively CDELT=1), then we'll rename the
    PC matrix keywords to CD. Because using the CD matrix is non-standard for
    Gnuastro, an explanation is added just after the WCS keywords, thoroughly
    explaining why we were forced to use the CD matrix.
    
    This fixes bug #58974.
---
 NEWS                         |   1 +
 doc/announce-acknowledge.txt |   1 +
 doc/gnuastro.texi            |   8 ++-
 lib/fits.c                   |  24 ++-----
 lib/gnuastro/wcs.h           |   3 +
 lib/wcs.c                    | 154 ++++++++++++++++++++++++++++++++++++++-----
 6 files changed, 153 insertions(+), 38 deletions(-)

diff --git a/NEWS b/NEWS
index 15517c6..c1608a4 100644
--- a/NEWS
+++ b/NEWS
@@ -121,6 +121,7 @@ See the end of the file for license conditions.
   bug #58835: Floating point errors when comparing pixel scale in Crop.
   bug #58898: Plain text string columns touching next, clear first character.
   bug #58901: Blank values for non-standard integer types in FITS tables.
+  bug #58974: WCS conversion not reasonable on processed TPV data.
 
 
 
diff --git a/doc/announce-acknowledge.txt b/doc/announce-acknowledge.txt
index 5633ca5..9f742ee 100644
--- a/doc/announce-acknowledge.txt
+++ b/doc/announce-acknowledge.txt
@@ -3,6 +3,7 @@ Alphabetically ordered list to acknowledge in the next release.
 Marjan Akbari
 Carlos Allende Prieto
 Leindert Boogaard
+Mark Calabretta
 Alexey Dokuchaev
 Raúl Infante Sainz
 Samane Raji
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 4263b56..b443683 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -22693,13 +22693,19 @@ number of coordinate representations found into the 
space that @code{nwcs}
 points to. Please see @code{gal_wcs_read_fitsptr} for more.
 @end deftypefun
 
-@deftypefun gal_wcs_write (struct wcsprm @code{*wcs}, char @code{*filename}, 
char @code{*extname}, gal_fits_list_key_t @code{*headers}, char 
@code{*program_string})
+@deftypefun void gal_wcs_write (struct wcsprm @code{*wcs}, char 
@code{*filename}, char @code{*extname}, gal_fits_list_key_t @code{*headers}, 
char @code{*program_string})
 Write the given WCS structure into the second extension of an empty FITS 
header.
 The first/primary extension will be empty like the default format of all 
Gnuastro outputs.
 When @code{extname!=NULL} it will be used as the FITS extension name.
 Any set of extra headers can also be written through the @code{headers} list 
and if @code{program_string!=NULL} it will be used in a commented keyword title 
just above the written version information.
 @end deftypefun
 
+@deftypefun void gal_wcs_write_in_fitsptr(fitsfile @code{*fptr}, struct wcsprm 
@code{*wcs})
+Convert the input @code{wcs} structure (keeping the WCS programmatically) into 
FITS keywords and write them into the given FITS file pointer.
+This is a relatively low-level function which assumes the FITS file has 
already been opened with CFITSIO.
+If you just want to write the WCS into an empty file, you can use 
@code{gal_wcs_write} (which internally calls this function after creating the 
FITS file and later closes it safely).
+@end deftypefun
+
 @deftypefun {struct wcsprm *} gal_wcs_copy (struct wcsprm @code{*wcs})
 Return a fully allocated (independent) copy of @code{wcs}.
 @end deftypefun
diff --git a/lib/fits.c b/lib/fits.c
index ce89c09..046d536 100644
--- a/lib/fits.c
+++ b/lib/fits.c
@@ -2034,12 +2034,12 @@ gal_fits_img_write_to_ptr(gal_data_t *input, char 
*filename)
 {
   void *blank;
   int64_t *i64;
+  char *u64key;
   fitsfile *fptr;
   uint64_t *u64, *u64f;
   long fpixel=1, *naxes;
-  char *wcsstr, *u64key;
   size_t i, ndim=input->ndim;
-  int nkeyrec, hasblank, status=0, datatype=0;
+  int hasblank, status=0, datatype=0;
   gal_data_t *i64data, *towrite, *block=gal_tile_block(input);
 
   /* Small sanity check. */
@@ -2174,25 +2174,11 @@ gal_fits_img_write_to_ptr(gal_data_t *input, char 
*filename)
   if(towrite->comment)
     fits_write_comment(fptr, towrite->comment, &status);
 
+
   /* If a WCS structure is present, write it in */
   if(towrite->wcs)
-    {
-      /* 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)
-        error(0, 0, "%s: WARNING: WCSLIB error, no WCS in output.\n"
-              "wcshdu ERROR %d: %s", __func__, status,
-              wcs_errmsg[status]);
-      else
-        gal_fits_key_write_wcsstr(fptr, wcsstr, nkeyrec);
-      status=0;
-    }
+    gal_wcs_write_in_fitsptr(fptr, towrite->wcs);
+
 
   /* Report any errors if we had any */
   free(naxes);
diff --git a/lib/gnuastro/wcs.h b/lib/gnuastro/wcs.h
index b76eabd..dede089 100644
--- a/lib/gnuastro/wcs.h
+++ b/lib/gnuastro/wcs.h
@@ -96,6 +96,9 @@ gal_wcs_write(struct wcsprm *wcs, char *filename,
               char *extname, gal_fits_list_key_t *headers,
               char *program_string);
 
+void
+gal_wcs_write_in_fitsptr(fitsfile *fptr, struct wcsprm *wcs);
+
 
 
 
diff --git a/lib/wcs.c b/lib/wcs.c
index 7123339..91ff028 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -50,6 +50,14 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 
 
+/* Static functions on for this file. */
+static void
+gal_wcs_to_cd(struct wcsprm *wcs);
+
+
+
+
+
 /*************************************************************
  ***********               Read WCS                ***********
  *************************************************************/
@@ -307,15 +315,76 @@ gal_wcs_read(char *filename, char *hdu, size_t hstartwcs,
  ***********               Write WCS               ***********
  *************************************************************/
 void
+gal_wcs_write_in_fitsptr(fitsfile *fptr, struct wcsprm *wcs)
+{
+  char *wcsstr;
+  int tpvdist, status=0, nkeyrec;
+
+  /* Prepare the main rotation matrix. Note that for TPV distortion, WCSLIB
+     versions 7.3 and before couldn't deal with the CDELT keys, so to be
+     safe, in such cases, we'll remove the effect of CDELT in the
+     'gal_wcs_to_cd' function. */
+  tpvdist=wcs->lin.disseq && !strcmp(wcs->lin.disseq->dtype[1], "TPV");
+  if( tpvdist ) gal_wcs_to_cd(wcs);
+  else          gal_wcs_decompose_pc_cdelt(wcs);
+
+  /* Clean up small errors in the PC matrix and CDELT values. */
+  gal_wcs_clean_errors(wcs);
+
+  /* Convert the WCS information to text. */
+  status=wcshdo(WCSHDO_safe, wcs, &nkeyrec, &wcsstr);
+  if(status)
+    error(0, 0, "%s: WARNING: WCSLIB error, no WCS in output.\n"
+          "wcshdu ERROR %d: %s", __func__, status,
+          wcs_errmsg[status]);
+  else
+    gal_fits_key_write_wcsstr(fptr, wcsstr, nkeyrec);
+  status=0;
+
+   /* WCSLIB is going to write PC+CDELT keywords in any case. But when we
+      have a TPV distortion, it is cleaner to use a CD matrix. Also,
+      including and before version 7.3, WCSLIB wouldn't convert coordinates
+      properly if the PC matrix is used with the TPV distortion. So to help
+      users with WCSLIB 7.3 or earlier, we need to conver the PC matrix to
+      CD. 'gal_wcs_to_cd' function already made sure that CDELT=1, so
+      effectively the CD matrix and PC matrix are equivalent, we just need
+      to convert the keyword names and delete the CDELT keywords. Note that
+      zero-valued PC/CD elements may not be present, so we'll manually set
+      'status' to zero and not let CFITSIO crash.*/
+  if(tpvdist)
+    {
+      status=0; fits_modify_name(fptr, "PC1_1", "CD1_1", &status);
+      status=0; fits_modify_name(fptr, "PC1_2", "CD1_2", &status);
+      status=0; fits_modify_name(fptr, "PC2_1", "CD2_1", &status);
+      status=0; fits_modify_name(fptr, "PC2_2", "CD2_2", &status);
+      status=0; fits_delete_str(fptr, "CDELT1", &status);
+      status=0; fits_delete_str(fptr, "CDELT2", &status);
+      status=0;
+      fits_write_comment(fptr, "The CD matrix is used instead of the "
+                         "PC+CDELT due to conflicts with TPV distortion "
+                         "in WCSLIB 7.3 (released on 2020/06/03) and "
+                         "ealier. By default Gnuastro will write "
+                         "PC+CDELT matrices because the rotation (PC) and "
+                         "pixel-scale (CDELT) are separate; providing "
+                         "more physically relevant metadata for human "
+                         "readers (PC+CDELT is also the default format "
+                         "of WCSLIB).", &status);
+    }
+}
+
+
+
+
+
+void
 gal_wcs_write(struct wcsprm *wcs, char *filename,
               char *extname, gal_fits_list_key_t *headers,
               char *program_string)
 {
-  char *wcsstr;
+  int status=0;
   size_t ndim=0;
   fitsfile *fptr;
   long *naxes=NULL;
-  int status=0, nkeyrec;
 
   /* Small sanity checks */
   if(wcs==NULL)
@@ -339,24 +408,12 @@ gal_wcs_write(struct wcsprm *wcs, char *filename,
   fits_delete_key(fptr, "COMMENT", &status);
   status=0;
 
+  /* If an extension name was requested, add it. */
   if(extname)
     fits_write_key(fptr, TSTRING, "EXTNAME", extname, "", &status);
 
-  /* Decompose the 'PCi_j' matrix and 'CDELTi' vector. */
-  gal_wcs_decompose_pc_cdelt(wcs);
-
-  /* Clean up small errors in the PC matrix and CDELT values. */
-  gal_wcs_clean_errors(wcs);
-
-  /* Convert the WCS information to text. */
-  status=wcshdo(WCSHDO_safe, wcs, &nkeyrec, &wcsstr);
-  if(status)
-    error(0, 0, "%s: WARNING: WCSLIB error, no WCS in output.\n"
-          "wcshdu ERROR %d: %s", __func__, status,
-          wcs_errmsg[status]);
-  else
-    gal_fits_key_write_wcsstr(fptr, wcsstr, nkeyrec);
-  status=0;
+  /* Write the WCS structure. */
+  gal_wcs_write_in_fitsptr(fptr, wcs);
 
   /* Write all the headers and the version information. */
   gal_fits_key_write_version_in_ptr(&headers, program_string, fptr);
@@ -1030,7 +1087,9 @@ gal_wcs_decompose_pc_cdelt(struct wcsprm *wcs)
   /* If there is on WCS, then don't do anything. */
   if(wcs==NULL) return;
 
-  /* The correction is only needed when the PC matrix is filled. */
+  /* The correction is only needed when the PC matrix is filled. Note that
+     internally, WCSLIB always uses the PC matrix, even when the input has
+     a CD matrix.*/
   if(wcs->pc)
     {
       /* Get the pixel scale. */
@@ -1072,6 +1131,65 @@ gal_wcs_decompose_pc_cdelt(struct wcsprm *wcs)
 
 
 
+/* Set the WCS structure to use the CD matrix. */
+static void
+gal_wcs_to_cd(struct wcsprm *wcs)
+{
+  size_t i, j, naxis;
+
+  /* If there is on WCS, then don't do anything. */
+  if(wcs==NULL) return;
+
+  /* 'wcs->altlin' identifies which rotation element is being used (PCi_j,
+     CDi_J or CROTAi). For PCi_j, the first bit will be 1 (==1), for CDi_j,
+     the second bit is 1 (==2) and for CROTAi, the third bit is 1 (==4). */
+  naxis=wcs->naxis;
+  switch(wcs->altlin)
+    {
+   /* PCi_j: Convert it to CDi_j. */
+    case 1:
+
+      /* Fill in the CD matrix and correct the PC and CDELT arrays. We have
+         to do this because ultimately, WCSLIB will be writing the PC and
+         CDELT keywords, even when 'altlin' is 2. So effectively we have to
+         multiply the PC and CDELT matrices, then set cdelt=1 in all
+         dimensions. This is actually how WCSLIB reads a FITS header with
+         only a CD matrix. */
+      for(i=0;i<naxis;++i)
+        {
+          for(j=0;j<naxis;++j)
+            wcs->cd[i*naxis+j] = wcs->pc[i*naxis+j] *= wcs->cdelt[i];
+          wcs->cdelt[i]=1;
+        }
+
+      /* Set the altlin to be the CD matrix and free the PC matrix. */
+      wcs->altlin=2;
+      break;
+
+    /* CDi_j: No need to do any conversion. */
+    case 2: return; break;
+
+    /* CROTAi: not yet supported. */
+    case 4:
+      error(0, 0, "%s: WARNING: Conversion of 'CROTAi' keywords to the CD "
+            "matrix is not yet supported (for lack of time!), please "
+            "contact us at %s to add this feature. But this may not cause a "
+            "problem at all, so please check if the output's WCS is "
+            "reasonable", __func__, PACKAGE_BUGREPORT);
+      break;
+
+    /* The value isn't supported! */
+    default:
+      error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix the "
+            "problem. The value %d for wcs->altlin isn't recognized",
+            __func__, PACKAGE_BUGREPORT, wcs->altlin);
+    }
+}
+
+
+
+
+
 /* The distance (along a great circle) on a sphere between two points
    is calculated here. Since the pixel sides are usually very small,
    we won't be using the direct formula: