[gnuastro-commits] master 9067a69: Library (wcs.h): access to WCSLIB's wcshdo() only within library

[Mohammad Akhlaghi]

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

    Library (wcs.h): access to WCSLIB's wcshdo() only within library
    
    Until now, the Crop and MakeProfiles programs directly called the 'wcshdo'
    function of WCSLIB (that is in charge of writing the FITS
    keywords). However, in time we have found some peculiarities with the WCS
    keyword writing, for example if there is a TPV distortion, the CD matrix
    should be used or some cleaning of the written parameters (sometimes
    happening due to floating point errors).
    
    With this commit a new function has been added to the 'wcs.h' library of
    Gnuastro: 'gal_wcs_write_wcsstr' and the only place that we call WCSLIB's
    'wcshdo' within the whole of Gnuastro is within this function. This greatly
    helps in modularity and unifying the output behavior of all programs.
    
    In the same spirit, the final CD/PC matrix correction (after the keywords
    were written into the FITS file) have been moved to
    'gal_fits_key_write_wcsstr' and the check for seeing if a CD matrix should
    be used has been moved into a function and used in the
    'gal_wcs_read_fitsptr' function (which is the only place we read WCS in
    Gnuastro).
    
    The necessity for all this was found after noting that the 'wcstoimg'
    operator in column arithmetic is giving a WCSLIB error when the image has
    TPV distortion (bug #60909). After contacting Mark Calabretta, he kindly
    noted that the problem is in not using the CD matrix with TPV. I then
    discovered that the tests we had added for this were not actually being
    used by the Crop or MakeProfiles programs. So instead of repeating those
    checks in these two programs (which would certainly cause many headaches in
    the future), I made the generic changes above and now we have a single
    place for reading/writing the WCS that will account for all existing/future
    corrections.
    
    This fixes bug #60909.
---
 NEWS                         |  12 +++--
 bin/crop/onecrop.c           |  10 +---
 bin/crop/ui.c                |   8 +---
 bin/mkprof/main.h            |   2 +-
 bin/mkprof/mkprof.c          |   2 +-
 bin/mkprof/ui.c              |  39 +++++++--------
 doc/announce-acknowledge.txt |   1 +
 doc/gnuastro.texi            |   7 +++
 lib/fits.c                   |  34 +++++++++++++
 lib/gnuastro/wcs.h           |   3 ++
 lib/wcs.c                    | 111 ++++++++++++++++++++++---------------------
 11 files changed, 132 insertions(+), 97 deletions(-)

diff --git a/NEWS b/NEWS
index a1475f0..c0998da 100644
--- a/NEWS
+++ b/NEWS
@@ -37,6 +37,8 @@ See the end of the file for license conditions.
   Library:
    - Arithmetic macros:
      - GAL_ARITHMETIC_OP_BOX_AROUND_ELLIPSE
+   - gal_wcs_write_wcsstr: wrapper of wcshdo of WCSLIB to simplify generic
+     access to that low-level operation (has some extra sanity checks).
 
 ** Removed features
 
@@ -56,15 +58,17 @@ See the end of the file for license conditions.
 ** Bugs fixed
   bug #60725: MakeCatalog doesn't put comment on --halfsumsb column.
   bug #60776: Radial profile script not using standard deviation image,
-              this bug was fixed by Zahra Sharbaf.
+              fixed by Zahra Sharbaf.
   bug #60778: Brightness error not NaN when all STD pixels are blank,
-              this bug was reported by Zahra Sharbaf.
+              reported by Zahra Sharbaf.
   bug #60826: Arithmetic won't delete existing file with tofile operators.
   bug #60881: Query segmentation fault when NED is called without a dataset.
   bug #60901: sort-by-night crash during make check due to stdintimeout,
-              this bug was reported by Zahra Sharbaf.
+              reported by Zahra Sharbaf.
   bug #60903: Increasing value of stdintimeout not effective beyond 1000000,
-              this bug was reported by Zahra Sharbaf.
+              reported by Zahra Sharbaf.
+  bug #60909: WCS coordinate conversion not working with TPV distortion,
+              fixed with the help of Mark Calabretta.
 
 
 
diff --git a/bin/crop/onecrop.c b/bin/crop/onecrop.c
index 470e966..f1c8bbd 100644
--- a/bin/crop/onecrop.c
+++ b/bin/crop/onecrop.c
@@ -613,15 +613,7 @@ onecrop_make_array(struct onecropparams *crp, long 
*fpixel_i,
   if(img->wcs)
     {
       /* Write the WCS title and common WCS information. */
-      if(fits_write_record(ofp, blankrec, &status))
-        gal_fits_io_error(status, NULL);
-      sprintf(titlerec, "%sWCS information", GAL_FITS_KEY_TITLE_START);
-      for(i=strlen(titlerec);i<79;++i)
-        titlerec[i]=' ';
-      fits_write_record(ofp, titlerec, &status);
-      for(i=0;i<img->nwcskeys-1;++i)
-        fits_write_record(ofp, &img->wcstxt[i*80], &status);
-      gal_fits_io_error(status, NULL);
+      gal_fits_key_write_wcsstr(ofp, img->wcs, img->wcstxt, img->nwcskeys);
 
       /* Correct the CRPIX keywords. */
       for(i=0;i<ndim;++i)
diff --git a/bin/crop/ui.c b/bin/crop/ui.c
index b22e4c9..892e101 100644
--- a/bin/crop/ui.c
+++ b/bin/crop/ui.c
@@ -1057,13 +1057,7 @@ ui_preparations(struct cropparams *p)
       img->wcs=gal_wcs_read_fitsptr(tmpfits, p->cp.wcslinearmatrix,
                                     p->hstartwcs, p->hendwcs, &img->nwcs);
       if(img->wcs)
-        {
-          gal_wcs_decompose_pc_cdelt(img->wcs);
-          status=wcshdo(0, img->wcs, &img->nwcskeys, &img->wcstxt);
-          if(status)
-            error(EXIT_FAILURE, 0, "wcshdo ERROR %d: %s", status,
-                  wcs_errmsg[status]);
-        }
+        img->wcstxt=gal_wcs_write_wcsstr(img->wcs, &img->nwcskeys);
       else
         if(p->mode==IMGCROP_MODE_WCS)
           error(EXIT_FAILURE, 0, "%s (hdu %s): the WCS structure is "
diff --git a/bin/mkprof/main.h b/bin/mkprof/main.h
index 7440f8c..d794946 100644
--- a/bin/mkprof/main.h
+++ b/bin/mkprof/main.h
@@ -188,7 +188,7 @@ struct mkprofparams
   pthread_cond_t     qready;  /* bq is ready to be written.               */
   pthread_mutex_t     qlock;  /* Mutex lock to change builtq.             */
   double          halfpixel;  /* Half pixel in oversampled image.         */
-  char           *wcsheader;  /* The WCS header information for main img. */
+  char              *wcsstr;  /* The WCS keywords derived from main img.  */
   int            wcsnkeyrec;  /* The number of keywords in the WCS header.*/
   char       *mergedimgname;  /* Name of merged image.                    */
   gal_data_t        *custom;  /* Table containing custom values.          */
diff --git a/bin/mkprof/mkprof.c b/bin/mkprof/mkprof.c
index 13b8cde..afcfaf1 100644
--- a/bin/mkprof/mkprof.c
+++ b/bin/mkprof/mkprof.c
@@ -162,7 +162,7 @@ saveindividual(struct mkonthread *mkp)
         crpix[i] = ((double *)(p->crpix->array))[i] - os*(mkp->fpixel_i[i]-1);
 
       /* Write the image. */
-      gal_fits_img_write_corr_wcs_str(ibq->image, filename, p->wcsheader,
+      gal_fits_img_write_corr_wcs_str(ibq->image, filename, p->wcsstr,
                                       p->wcsnkeyrec, crpix, NULL,
                                       PROGRAM_NAME);
     }
diff --git a/bin/mkprof/ui.c b/bin/mkprof/ui.c
index d12ede4..acbc2d0 100644
--- a/bin/mkprof/ui.c
+++ b/bin/mkprof/ui.c
@@ -1374,9 +1374,9 @@ static void
 ui_prepare_canvas(struct mkprofparams *p)
 {
   float *f, *ff;
+  int setshift=0;
   long width[3]={1,1,1};
   size_t tndim, *tdsize;
-  int status=0, setshift=0;
   double truncr, semiaxes[3], euler_deg[3];
   size_t i, nshift=0, *dsize=NULL, ndim_counter;
 
@@ -1528,19 +1528,6 @@ ui_prepare_canvas(struct mkprofparams *p)
       if(p->out->unit==NULL)
         gal_checkset_allocate_copy("Brightness", &p->out->unit);
     }
-
-
-  /* When individual mode is requested, write the WCS structure to a header
-     string to speed up the process: if we don't do it here, this process
-     will be necessary on every individual profile's output. So it is much
-     more efficient done once here. */
-  if(p->individual && p->wcs)
-    {
-      status=wcshdo(WCSHDO_safe, p->wcs, &p->wcsnkeyrec, &p->wcsheader);
-      if(status)
-        error(EXIT_FAILURE, 0, "wcshdo error %d: %s", status,
-              wcs_errmsg[status]);
-    }
 }
 
 
@@ -1859,11 +1846,22 @@ ui_preparations(struct mkprofparams *p)
   else
     ui_prepare_canvas(p);
 
-  /* Read the (possible) RA/Dec inputs into X and Y for the builder.  NOTE:
-     It may happen that there are no input columns, in that case, just
-     ignore this step.*/
-  if(p->wcs && p->num)
-    ui_finalize_coordinates(p);
+  /* Preparations now that we have WCS (if any was given in any way: either
+     from a background or from options).  */
+  if(p->wcs)
+    {
+      /* Read the (possible) RA/Dec inputs into X and Y for the builder.
+         NOTE: It may happen that there are no input columns, in that case,
+         just ignore this step.*/
+      if(p->num)
+        ui_finalize_coordinates(p);
+
+      /* If individual mode is activated, write the WCS as a string here
+         (earlier than needed). This is because it will be necessary for
+         every individual profile, but it will identical (except for the
+         'CRPIX's that will be changed). */
+      p->wcsstr=gal_wcs_write_wcsstr(p->wcs, &p->wcsnkeyrec);
+    }
 
   /* Prepare the random number generator. */
   p->rng=gal_checkset_gsl_rng(p->envseed, &p->rng_name, &p->rng_seed);
@@ -2064,8 +2062,7 @@ ui_free_report(struct mkprofparams *p, struct timeval *t1)
     }
 
   /* Free the WCS headers string that was defined for individual mode. */
-  if(p->individual)
-    free(p->wcsheader);
+  if(p->wcsstr) free(p->wcsstr);
 
   /* Free the random number generator: */
   gsl_rng_free(p->rng);
diff --git a/doc/announce-acknowledge.txt b/doc/announce-acknowledge.txt
index bf1d900..90c6cdf 100644
--- a/doc/announce-acknowledge.txt
+++ b/doc/announce-acknowledge.txt
@@ -1,5 +1,6 @@
 Alphabetically ordered list to acknowledge in the next release.
 
+Mark Calabretta
 Leslie Hunt
 Juan Molina Tobar
 Zahra Sharbaf
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 12e0678..157c491 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -26810,6 +26810,13 @@ For example if @code{CTYPE1} is @code{RA---TAN}, the 
string that function return
 Recall that the second component of @code{CTYPEi} contains the type of 
projection.
 @end deftypefun
 
+@deftypefun {char *} gal_wcs_write_wcsstr (struct wcsprm @code{*wcs}, int 
@code{*nkeyrec})
+Return an allocated string which contains the respective FITS keywords for the 
given WCS structure into it.
+The number of keywords is written in the space pointed by @code{nkeyrec}.
+Each FITS keyword is 80 characters wide (according to the FITS standard), and 
the next one is placed immediately after it, so the full string has 
@code{80*nkeyrec} bytes.
+The output of this function can later be written into an opened FITS file 
using @code{gal_fits_key_write_wcsstr} (see @ref{FITS header keywords}).
+@end deftypefun
+
 @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.
diff --git a/lib/fits.c b/lib/fits.c
index 140cf09..839671a 100644
--- a/lib/fits.c
+++ b/lib/fits.c
@@ -1895,6 +1895,9 @@ gal_fits_key_write_wcsstr(fitsfile *fptr, struct wcsprm 
*wcs,
   int status=0;
   char *keystart;
 
+  /* If the 'wcsstr' string is empty, don't do anything, simply return. */
+  if(wcsstr==NULL) return;
+
   /* Write the title. */
   gal_fits_key_write_title_in_ptr("World Coordinate System (WCS)", fptr);
 
@@ -1917,6 +1920,37 @@ gal_fits_key_write_wcsstr(fitsfile *fptr, struct wcsprm 
*wcs,
         }
     }
   gal_fits_io_error(status, NULL);
+
+   /* WCSLIB is going to write PC+CDELT keywords in any case. But when we
+      have a TPV, TNX or ZPX distortion, it is cleaner to use a CD matrix
+      (WCSLIB can't convert coordinates properly if the PC matrix is used
+      with the TPV distortion). So to help users avoid the potential
+      problems with WCSLIB, upon reading the WCS structure, in such cases,
+      we set CDELTi=1.0 and 'altlin=2' (signifying that the CD matrix
+      should be used). Therefore, 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 to avoid CFITSIO from
+      crashing. */
+  if(wcs && wcs->altlin==2)
+    {
+      status=0; fits_delete_str(fptr, "CDELT1", &status);
+      status=0; fits_delete_str(fptr, "CDELT2", &status);
+      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);
+      if(wcs->naxis==3)
+        {
+          status=0; fits_delete_str(fptr, "CDELT3", &status);
+          status=0; fits_modify_name(fptr, "PC1_3", "CD1_3", &status);
+          status=0; fits_modify_name(fptr, "PC2_3", "CD2_3", &status);
+          status=0; fits_modify_name(fptr, "PC3_1", "CD3_1", &status);
+          status=0; fits_modify_name(fptr, "PC3_2", "CD3_2", &status);
+          status=0; fits_modify_name(fptr, "PC3_3", "CD3_3", &status);
+        }
+      status=0;
+    }
 }
 
 
diff --git a/lib/gnuastro/wcs.h b/lib/gnuastro/wcs.h
index dfe430a..4abe42b 100644
--- a/lib/gnuastro/wcs.h
+++ b/lib/gnuastro/wcs.h
@@ -125,6 +125,9 @@ gal_wcs_write(struct wcsprm *wcs, char *filename,
               char *extname, gal_fits_list_key_t *headers,
               char *program_string);
 
+char *
+gal_wcs_write_wcsstr(struct wcsprm *wcs, int *nkeyrec);
+
 void
 gal_wcs_write_in_fitsptr(fitsfile *fptr, struct wcsprm *wcs);
 
diff --git a/lib/wcs.c b/lib/wcs.c
index b29db3d..7a64ca5 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -124,6 +124,24 @@ wcs_read_correct_pc_cd(struct wcsprm *wcs)
 
 
 
+/* For the TPV, TNX and ZPX distortions, WCSLIB can't deal with the CDELT
+   keys properly and its better to use the CD matrix instead, this function
+   will check for this and return 1 if a CD matrix should be used
+   (over-riding the user's desired matrix if necessary). */
+static int
+wcs_use_cd_for_distortion(struct wcsprm *wcs)
+{
+  return ( wcs
+           && wcs->lin.disseq
+           && ( !strcmp(   wcs->lin.disseq->dtype[1], "TPV")
+                || !strcmp(wcs->lin.disseq->dtype[1], "TNX")
+                || !strcmp(wcs->lin.disseq->dtype[1], "ZPX") ) );
+}
+
+
+
+
+
 /* Read the WCS information from the header. Unfortunately, WCS lib is
    not thread safe, so it needs a mutex. In case you are not using
    multiple threads, just pass a NULL pointer as the mutex.
@@ -369,10 +387,12 @@ gal_wcs_read_fitsptr(fitsfile *fptr, int linearmatrix, 
size_t hstartwcs,
         }
     }
 
-  /* If the user wants a CD linear matrix, do the conversion here,
-     otherwise, make sure the PC matrix is used. */
-  if(linearmatrix==GAL_WCS_LINEAR_MATRIX_CD) gal_wcs_to_cd(wcs);
-  else                          gal_wcs_decompose_pc_cdelt(wcs);
+  /* If the distortion requires a CD matrix, or if the user wants it, do
+     the conversion here, otherwise, make sure the PC matrix is used. */
+  if(wcs_use_cd_for_distortion(wcs) \
+     || linearmatrix==GAL_WCS_LINEAR_MATRIX_CD)
+    gal_wcs_to_cd(wcs);
+  else gal_wcs_decompose_pc_cdelt(wcs);
 
   /* Clean up and return. */
   status=0;
@@ -512,70 +532,53 @@ gal_wcs_dimension_name(struct wcsprm *wcs, size_t 
dimension)
 /*************************************************************
  ***********               Write WCS               ***********
  *************************************************************/
-void
-gal_wcs_write_in_fitsptr(fitsfile *fptr, struct wcsprm *wcs)
+char *
+gal_wcs_write_wcsstr(struct wcsprm *wcs, int *nkeyrec)
 {
   char *wcsstr;
-  int cdfordist, status=0, nkeyrec;
-
-  /* For the TPV, TNX and ZPX distortions, WCSLIB can't deal with the CDELT
-     keys properly and its better to use the CD matrix instead, so we'll
-     use the 'gal_wcs_to_cd' function. */
-  cdfordist = ( wcs->lin.disseq
-                && ( !strcmp(   wcs->lin.disseq->dtype[1], "TPV")
-                     || !strcmp(wcs->lin.disseq->dtype[1], "TNX")
-                     || !strcmp(wcs->lin.disseq->dtype[1], "ZPX") ) );
+  int status=0;
 
   /* Finalize the linear transformation matrix. Note that some programs may
      have worked on the WCS. So even if 'altlin' is already 2, we'll just
      ensure that the final matrix is CD here. */
-  if(wcs->altlin==2 || cdfordist) gal_wcs_to_cd(wcs);
-  else                            gal_wcs_decompose_pc_cdelt(wcs);
+  if(wcs->altlin==2 || wcs_use_cd_for_distortion(wcs)) 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);
+  /* Convert the WCS information to text. If there was an error, then free
+     any allocated space and put zero on 'nkeyrec'. */
+  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, wcs, wcsstr, nkeyrec);
-      free(wcsstr);
+      error(0, 0, "%s: WARNING: WCSLIB error, no WCS in output.\n"
+            "wcshdo ERROR %d: %s", __func__, status, wcs_errmsg[status]);
+      if(wcsstr) free(wcsstr);
+      *nkeyrec=0;
+      wcsstr=NULL;
     }
-  status=0;
 
-   /* WCSLIB is going to write PC+CDELT keywords in any case. But when we
-      have a TPV, TNX or ZPX distortion, it is cleaner to use a CD matrix
-      (WCSLIB can't convert coordinates properly if the PC matrix is used
-      with the TPV distortion). So to help users avoid the potential
-      problems with WCSLIB. 'gal_wcs_to_cd' function already made sure that
-      CDELTi=1.0, 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 to avoid CFITSIO from
-      crashing. */
-  if(wcs->altlin==2)
-    {
-      status=0; fits_delete_str(fptr, "CDELT1", &status);
-      status=0; fits_delete_str(fptr, "CDELT2", &status);
-      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);
-      if(wcs->naxis==3)
-        {
-          status=0; fits_delete_str(fptr, "CDELT3", &status);
-          status=0; fits_modify_name(fptr, "PC1_3", "CD1_3", &status);
-          status=0; fits_modify_name(fptr, "PC2_3", "CD2_3", &status);
-          status=0; fits_modify_name(fptr, "PC3_1", "CD3_1", &status);
-          status=0; fits_modify_name(fptr, "PC3_2", "CD3_2", &status);
-          status=0; fits_modify_name(fptr, "PC3_3", "CD3_3", &status);
-        }
-      status=0;
-    }
+  /* Return the string. */
+  return wcsstr;
+}
+
+
+
+
+
+void
+gal_wcs_write_in_fitsptr(fitsfile *fptr, struct wcsprm *wcs)
+{
+  char *wcsstr;
+  int nkeyrec=0;
+
+  /* Convert the WCS structure into FITS keywords as a string. */
+  wcsstr=gal_wcs_write_wcsstr(wcs, &nkeyrec);
+
+  /* Write the keywords into the FITS file. */
+  gal_fits_key_write_wcsstr(fptr, wcs, wcsstr, nkeyrec);
+  free(wcsstr);
 }