[gnuastro-commits] master 808c95d: Fits: default size with --wcsdistortion=SIP and no input data

[Mohammad Akhlaghi]

branch: master
commit 808c95dc56baf023928eeab3edf8bc6e3f572de0
Author: Sachin Kumar Singh <sachinkumarsingh092@gmail.com>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>

    Fits: default size with --wcsdistortion=SIP and no input data
    
    Until now, when the Fits program was called with '--wcsdistortion=SIP'
    (from a TPV-based WCS) and the input HDU didn't have any data, it would
    crash, complaining that a default size is necessary.
    
    With this commit, a set of steps have been added in the function called by
    '--wcsdistortion' to allow a default size when the input HDU has no
    size.
    
    In the process, a bug was fixed in reading the keywords by
    `wcsdistortion_get_tpvparams`. The problem was fixed by reading directly
    from 'disseq' parameter in wcs, which also made its reading style similar
    to that of `wcsdistortion_get_sipparams`.
---
 bin/fits/keywords.c | 38 +++++++++++++++++++++-----
 doc/gnuastro.texi   |  4 ++-
 lib/gnuastro/wcs.h  |  3 ++-
 lib/wcs.c           |  6 ++++-
 lib/wcsdistortion.c | 78 ++++++++++++++++++++++++++++++-----------------------
 5 files changed, 87 insertions(+), 42 deletions(-)

diff --git a/bin/fits/keywords.c b/bin/fits/keywords.c
index a4e47f3..c941539 100644
--- a/bin/fits/keywords.c
+++ b/bin/fits/keywords.c
@@ -442,9 +442,10 @@ keywords_distortion_wcs(struct fitsparams *p)
 {
   int nwcs;
   size_t ndim, *insize;
-  gal_data_t *data=NULL;
   char *suffix, *output;
+  gal_data_t *data=NULL;
   struct wcsprm *inwcs, *outwcs;
+  size_t *dsize, defaultsize[2]={2000,2000};
 
   /* If the extension has any data, read it, otherwise just make an empty
      array. */
@@ -461,11 +462,36 @@ keywords_distortion_wcs(struct fitsparams *p)
                                p->cp.quietmmap);
     }
 
-  /* Read the input WCS structure and convert it to the desired output
-     distortion. */
+  /* Read the input's WCS and make sure one exists. */
   inwcs=gal_wcs_read(p->filename, p->cp.hdu, 0, 0, &nwcs);
-  outwcs=gal_wcs_distortion_convert(inwcs, p->distortionid,
-                                    data?data->dsize:NULL);
+  if(inwcs==NULL)
+    error(EXIT_FAILURE, 0, "%s (hdu %s): doesn't have any WCS structure "
+          "for converting its distortion",
+          p->filename, p->cp.hdu);
+
+  /* In case there is no dataset and the conversion is between TPV to SIP,
+     we need to set a default size and use that for the conversion, but we
+     also need to warn the user. */
+  if(data==NULL)
+    {
+      if( !p->cp.quiet
+          && gal_wcs_distortion_identify(inwcs)==GAL_WCS_DISTORTION_TPV
+          && p->distortionid==GAL_WCS_DISTORTION_SIP )
+        error(0, 0, "no data associated with WCS for distortion "
+              "converstion.\n\n"
+              "The requested conversion can't be done analytically, so a "
+              "solution has to be found by fitting the parameters over a "
+              "grid of pixels. We will use a default grid of %zux%zu pixels "
+              "and will proceed with the conversion. But it would be more "
+              "accurate if it is the size of the image that this WCS is "
+              "associated with",
+              defaultsize[1], defaultsize[0]);
+      dsize=defaultsize;
+    }
+  else dsize=data->dsize;
+
+  /* Do the conversion. */
+  outwcs=gal_wcs_distortion_convert(inwcs, p->distortionid, dsize);
 
   /* Set the output filename. */
   if(p->cp.output)
@@ -491,7 +517,7 @@ keywords_distortion_wcs(struct fitsparams *p)
       gal_data_free(data);
     }
   else
-    gal_wcs_write(outwcs, output, NULL, PROGRAM_NAME);
+    gal_wcs_write(outwcs, output, p->wcsdistortion, NULL, PROGRAM_NAME);
 
   /* Clean up. */
   wcsfree(inwcs);
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 4b22347..0c17cf6 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -22192,9 +22192,11 @@ 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}, 
gal_fits_list_key_t @code{*headers}, char @code{*program_string})
+@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})
 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 {struct wcsprm *} gal_wcs_copy (struct wcsprm @code{*wcs})
diff --git a/lib/gnuastro/wcs.h b/lib/gnuastro/wcs.h
index 5821e3c..b76eabd 100644
--- a/lib/gnuastro/wcs.h
+++ b/lib/gnuastro/wcs.h
@@ -93,7 +93,8 @@ gal_wcs_read(char *filename, char *hdu, size_t hstartwcs,
  *************************************************************/
 void
 gal_wcs_write(struct wcsprm *wcs, char *filename,
-              gal_fits_list_key_t *headers, char *program_string);
+              char *extname, gal_fits_list_key_t *headers,
+              char *program_string);
 
 
 
diff --git a/lib/wcs.c b/lib/wcs.c
index a19d00d..1170b01 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -276,7 +276,8 @@ gal_wcs_read(char *filename, char *hdu, size_t hstartwcs,
  *************************************************************/
 void
 gal_wcs_write(struct wcsprm *wcs, char *filename,
-              gal_fits_list_key_t *headers, char *program_string)
+              char *extname, gal_fits_list_key_t *headers,
+              char *program_string)
 {
   char *wcsstr;
   size_t ndim=0;
@@ -306,6 +307,9 @@ gal_wcs_write(struct wcsprm *wcs, char *filename,
   fits_delete_key(fptr, "COMMENT", &status);
   status=0;
 
+  if(extname)
+    fits_write_key(fptr, TSTRING, "EXTNAME", extname, "", &status);
+
   /* Decompose the 'PCi_j' matrix and 'CDELTi' vector. */
   gal_wcs_decompose_pc_cdelt(wcs);
 
diff --git a/lib/wcsdistortion.c b/lib/wcsdistortion.c
index 56a3b5a..31212f7 100644
--- a/lib/wcsdistortion.c
+++ b/lib/wcsdistortion.c
@@ -74,59 +74,72 @@ static void
 wcsdistortion_get_tpvparams(struct wcsprm *wcs, double cd[2][2],
                             double *pv1, double *pv2)
 {
-  size_t pv_m=0;
-  size_t i, j, k,index=0;
+  const char *cp;
+  double *temp_cd=NULL;
+  struct dpkey  *keyp=NULL;
+  size_t i, j, k=0, index=0;
+  struct disprm *disseq=NULL;
 
   /* Make sure a WCS structure is actually given. */
   if(wcs==NULL)
     error(EXIT_FAILURE, 0, "%s: input WCS structure is NULL", __func__);
 
-  for(i=0,k=0; i<2; ++i)
+  /* For easy reading. */
+  disseq=wcs->lin.disseq;
+  keyp=disseq->dp;
+
+
+  /* Fill the 2-times allocated CD array (cd[][]). Note that the required
+    CD matix is extracted using the `gal_wcs_wrap_matrix` as a single
+    allocated array (temp_cd[]), that is then used to fill cd[][]. */
+  temp_cd=gal_wcs_warp_matrix(wcs);
+  for(i=0; i<2; ++i)
     for(j=0; j<2; ++j)
       {
         /* If a value is present store it, else store 0.*/
-        if((wcs->cd)[i] != 0)   cd[i][j]=(wcs->cd)[k++];
-        else                    cd[i][j]=0;
+        if(temp_cd[k] != 0)   cd[i][j]=temp_cd[k++];
+        else                  cd[i][j]=0;
 
         /* For a check:
-        printf("CD%ld_%ld\t%.10lf\n", i+1, j+1, cd[i][j]);
+        printf("CD%ld_%ld\t%.10E\n", i+1, j+1, cd[i][j]);
         */
       }
 
-  for(j=0; j < wcs->npvmax; ++j)
+  for(i=0; i<disseq->ndp; ++i, ++keyp)
     {
-      if(wcs->pv[j].i == 1)
+      /* For axis1. */
+      if (keyp->j == 1)
         {
-          /*pv_m is used to check the index in the header.*/
-          pv_m=wcs->pv[j].m;
-
-          /* `index` is the index of the pv* array.*/
-          index = pv_m;
+          /* Ignore any missing keyvalues. */
+          if ( keyp->field == NULL ) continue;
+          cp = strchr(keyp->field, '.') + 1;
+          if (strncmp(cp, "TPV.", 4) != 0) continue;
+          cp += 4;
 
-          if( wcs->pv[pv_m].value != 0 && index == pv_m )
-            pv1[index]=wcs->pv[j].value;
-          else
-            pv1[index]=0;
+          sscanf(cp, "%ld", &index);
+          pv1[index]=disseq->dp[i].value.f;
 
           /* For a check
-          printf("PV1_%d\t%.10f\n", pv_m, pv1[k]);
+          printf("PV1_%ld\t%.10f\n", index, pv1[index]);
           */
         }
-      else if(wcs->pv[j].i == 2)
+
+      /* For axis2. */
+      else if (keyp->j == 2)
         {
-          /*pv_m is used to check the index in the header.*/
-          pv_m=wcs->pv[j].m;
+          /* Ignore any missing keyvalues. */
+          if ( keyp->field == NULL ) continue;
+
+          cp = strchr(keyp->field, '.') + 1;
+          if (strncmp(cp, "TPV.", 4) != 0) continue;
+          cp += 4;
 
-          /* `index` is the index of the pv* array.*/
-          index = pv_m;
+          sscanf(cp, "%ld", &index);
 
-          if( wcs->pv[pv_m].value != 0 && index == pv_m )
-            pv2[index]=wcs->pv[j].value;
-          else
-            pv2[index]=0;
+          pv2[index]=disseq->dp[i].value.f;
 
           /* For a check
-          printf("PV2_%d\t%.10f\n", pv_m, pv2[k]);
+          printf("PV2_%ld\t%.10f\n", index, pv2[index]);
           */
         }
       else
@@ -513,7 +526,7 @@ wcsdistortion_intermidate_tpveq(double cd[2][2], double 
*pv1, double *pv2,
 
   /* For a check:
   {
-    size i, j;
+    size_t i, j;
     for(i=0; i<=4; ++i)
       for(j=0;j<=4;++j)
        {
@@ -1063,14 +1076,13 @@ wcsdistortion_calcsip(size_t axis, size_t m, size_t n, 
double tpvu[8][8],
   double sip_coeff;
   if(     axis == 1) sip_coeff=tpvu[m][n];
   else if(axis == 2) sip_coeff=tpvv[m][n];
-  else error(EXIT_FAILURE, 0, "%s: axis does not exists! ", __func__);
+  else error(EXIT_FAILURE, 0, "%s: axis %zu does not exists! ",
+             __func__, axis);
 
   if(      (axis == 1) && (m == 1) && (n == 0) ) sip_coeff = sip_coeff - 1.0;
   else if( (axis == 2) && (m == 0) && (n == 1) ) sip_coeff = sip_coeff - 1.0;
-  else error(EXIT_FAILURE, 0, "%s: axis does not exists! ", __func__);
 
   return sip_coeff;
-
 }
 
 
@@ -1432,7 +1444,7 @@ wcsdistortion_add_sipkeywords(struct wcsprm *wcs, size_t 
*fitsize,
   int size = wcs->naxis;
   size_t a_order=0, b_order=0;
   size_t ap_order=0, bp_order=0;
-  size_t m, n, num=0, numkey=100;
+  size_t m, n, num=0, numkey=200;
   double ap_coeff[5][5], bp_coeff[5][5];
   char *fullheader, fmt[50], sipkey[8], keyaxis[9], pcaxis[10];