[gnuastro-commits] master dc47d2a: Table: convert WCS and image coordinates in output columns

[Mohammad Akhlaghi]

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

    Table: convert WCS and image coordinates in output columns
    
    Until now, Table would just print the desired column values without
    actually changing their values.
    
    However, some common changing of values is often required during data
    analysis. With this commit, Table now has the ability to change its column
    values and print new ones. In particular, the two `wcstoimg' and `imgtowcs'
    have been defined to convert coordinates when creating the output.
---
 NEWS                    |   6 +++
 bin/table/args.h        |  26 ++++++++++
 bin/table/asttable.conf |   3 ++
 bin/table/main.h        |   6 +++
 bin/table/table.c       |  72 +++++++++++++++++++++++++++
 bin/table/ui.c          | 126 ++++++++++++++++++++++++++++++++++++++++++++++--
 bin/table/ui.h          |   6 ++-
 doc/gnuastro.texi       |  41 ++++++++++++++++
 8 files changed, 279 insertions(+), 7 deletions(-)

diff --git a/NEWS b/NEWS
index 1f656bd..c0cc5a9 100644
--- a/NEWS
+++ b/NEWS
@@ -27,6 +27,12 @@ See the end of the file for license conditions.
      the book under `--sigclip-median' for a nice use case.
 
   Table:
+   - As part of printing output columns, Table can now convert Image to WCS
+     coordinates and vice-versa. For example if the input catalog has
+     atleast an `ID' column and two `RA' and `DEC' columns the following
+     `-cID,RA,DEC -cwcstoimg(RA,DEC) --wcsfile=a.fits' options will print 5
+     columns, where the last two columns are the image coordinates based on
+     the WCS in `a.fits'.
    --head: Only output the given number of rows from the top of columns.
    --tail: Only output the given number of rows from the botoom of columns.
 
diff --git a/bin/table/args.h b/bin/table/args.h
index 76e5478..144616c 100644
--- a/bin/table/args.h
+++ b/bin/table/args.h
@@ -44,6 +44,32 @@ struct argp_option program_options[] =
       GAL_OPTIONS_NOT_MANDATORY,
       GAL_OPTIONS_NOT_SET
     },
+    {
+      "wcsfile",
+      UI_KEY_WCSFILE,
+      "STR",
+      0,
+      "File with WCS if conversion is requested.",
+      GAL_OPTIONS_GROUP_INPUT,
+      &p->wcsfile,
+      GAL_TYPE_STRING,
+      GAL_OPTIONS_RANGE_ANY,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
+    },
+    {
+      "wcshdu",
+      UI_KEY_WCSHDU,
+      "STR",
+      0,
+      "HDU in file with WCS for conversion.",
+      GAL_OPTIONS_GROUP_INPUT,
+      &p->wcshdu,
+      GAL_TYPE_STRING,
+      GAL_OPTIONS_RANGE_ANY,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
+    },
 
 
 
diff --git a/bin/table/asttable.conf b/bin/table/asttable.conf
index a7a1a97..4660480 100644
--- a/bin/table/asttable.conf
+++ b/bin/table/asttable.conf
@@ -18,3 +18,6 @@
 # permitted in any medium without royalty provided the copyright notice and
 # this notice are preserved.  This file is offered as-is, without any
 # warranty.
+
+# Inputs
+ wcshdu        1
\ No newline at end of file
diff --git a/bin/table/main.h b/bin/table/main.h
index 53d733e..005424b 100644
--- a/bin/table/main.h
+++ b/bin/table/main.h
@@ -53,6 +53,8 @@ struct tableparams
   /* From command-line */
   struct gal_options_common_params cp; /* Common parameters.            */
   char              *filename;  /* Input filename.                      */
+  char               *wcsfile;  /* File with WCS.                       */
+  char                *wcshdu;  /* HDU in file with WCS.                */
   gal_list_str_t     *columns;  /* List of given columns.               */
   uint8_t         information;  /* ==1: only print FITS information.    */
   uint8_t     colinfoinstdout;  /* ==1: print column metadata in CL.    */
@@ -64,12 +66,16 @@ struct tableparams
 
   /* Internal. */
   gal_data_t           *table;  /* Linked list of output table columns. */
+  struct wcsprm          *wcs;  /* WCS structure for conversion.        */
+  int                    nwcs;  /* Number of WCS structures.            */
   gal_data_t      *allcolinfo;  /* Information of all the columns.      */
   gal_data_t         *sortcol;  /* Column to define a sorting.          */
   struct list_range *rangecol;  /* Column to define a range.            */
   uint8_t            freesort;  /* If the sort column should be freed.  */
   uint8_t          *freerange;  /* If the range column should be freed. */
   uint8_t              sortin;  /* If the sort column is in the output. */
+  size_t             wcstoimg;  /* Output column no, for conversion.    */
+  size_t             imgtowcs;  /* Output column no, for conversion.    */
   time_t              rawtime;  /* Starting time of the program.        */
 };
 
diff --git a/bin/table/table.c b/bin/table/table.c
index 2f4e5e7..50b83c7 100644
--- a/bin/table/table.c
+++ b/bin/table/table.c
@@ -31,6 +31,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 #include <gsl/gsl_heapsort.h>
 
+#include <gnuastro/wcs.h>
 #include <gnuastro/fits.h>
 #include <gnuastro/table.h>
 #include <gnuastro/qsort.h>
@@ -321,6 +322,73 @@ table_head_tail(struct tableparams *p)
 
 
 
+
+
+static void
+table_unit_conversion(struct tableparams *p, int w0i1)
+{
+  int isfirstcol;
+  gal_data_t *t1, *t2, *tmp, *before, *after;
+  size_t i, startcol = w0i1 ? p->imgtowcs : p->wcstoimg;
+
+  /* Go to the column that we need to convert. */
+  i=0;
+  before=p->table;
+  for(tmp=p->table; tmp!=NULL; tmp=tmp->next)
+    {
+      if( i++ == startcol )
+        {
+          after=tmp->next->next;
+          break;
+        }
+      before=tmp;
+    }
+  isfirstcol=before==p->table;
+
+  /* Make sure the types are double-precision floating point. NOTE that
+     since we are freeing afterwards, the second column needs to be
+     read first.*/
+  t2=gal_data_copy_to_new_type_free(tmp->next, GAL_TYPE_FLOAT64);
+  t1=gal_data_copy_to_new_type_free(tmp,       GAL_TYPE_FLOAT64);
+
+  /* Do the conversion. */
+  t1->next=t2;
+  t2->next=NULL;
+  if(w0i1) gal_wcs_img_to_world(t1, p->wcs, 1);
+  else     gal_wcs_world_to_img(t1, p->wcs, 1);
+
+  /* In image mode, double-precision floating point is too much. */
+  if(w0i1==0)
+    {
+      t1=gal_data_copy_to_new_type_free(t1, GAL_TYPE_FLOAT32);
+      t2=gal_data_copy_to_new_type_free(t2, GAL_TYPE_FLOAT32);
+    }
+
+  /* Put them back into the output table. */
+  t1->next=t2;
+  t2->next=after;
+  if(isfirstcol) p->table=t1; else before->next=t1;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
 /**************************************************************/
 /***************       Top function         *******************/
 /**************************************************************/
@@ -337,6 +405,10 @@ table(struct tableparams *p)
   if(p->head!=GAL_BLANK_SIZE_T || p->tail!=GAL_BLANK_SIZE_T)
     table_head_tail(p);
 
+  /* If unit conversion is requested, do it. */
+  if(p->wcstoimg!=GAL_BLANK_SIZE_T) table_unit_conversion(p, 0);
+  if(p->imgtowcs!=GAL_BLANK_SIZE_T) table_unit_conversion(p, 1);
+
   /* Write the output. */
   gal_table_write(p->table, NULL, p->cp.tableformat, p->cp.output,
                   "TABLE", p->colinfoinstdout);
diff --git a/bin/table/ui.c b/bin/table/ui.c
index 8f7a032..392a94e 100644
--- a/bin/table/ui.c
+++ b/bin/table/ui.c
@@ -28,6 +28,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <stdio.h>
 #include <string.h>
 
+#include <gnuastro/wcs.h>
 #include <gnuastro/fits.h>
 #include <gnuastro/table.h>
 #include <gnuastro/pointer.h>
@@ -118,6 +119,8 @@ ui_initialize_options(struct tableparams *p,
   /* Program-specific initialization. */
   p->head                = GAL_BLANK_SIZE_T;
   p->tail                = GAL_BLANK_SIZE_T;
+  p->wcstoimg            = GAL_BLANK_SIZE_T;
+  p->imgtowcs            = GAL_BLANK_SIZE_T;
 
   /* Modify common options. */
   for(i=0; !gal_options_is_last(&cp->coptions[i]); ++i)
@@ -432,17 +435,93 @@ ui_print_info_exit(struct tableparams *p)
 
 
 
+/* WCS <-> Image conversion */
+static void
+ui_wcs_conversion(struct tableparams *p, char *instr,
+                  gal_list_str_t **wcstoimg_ptr,
+                  gal_list_str_t **imgtowcs_ptr,
+                  gal_list_str_t **new)
+{
+  const char delimiter[]="/";
+  char *c1, *c2, *c3, *optname, *saveptr;
+
+  /* Set the basic option properties. */
+  instr[8]='\0';
+  optname = instr;
+
+  /* If a WCS hasn't been read yet, read it.*/
+  if(p->wcs==NULL)
+    {
+      /* A small sanity check. */
+      if(p->wcsfile==NULL || p->wcshdu==NULL)
+        error(EXIT_FAILURE, 0, "`--wcsfile' and `--wcshdu' are necessary for "
+              "the `%s' conversion", optname);
+
+      /* Read the WCS. */
+      p->wcs=gal_wcs_read(p->wcsfile, p->wcshdu, 0, 0, &p->nwcs);
+      if(p->wcs==NULL)
+        error(EXIT_FAILURE, 0, "%s (hdu: %s): no WCS could be read by "
+              "WCSLIB", p->wcsfile, p->wcshdu);
+    }
+
+  /* Make sure this conversion is only requested once. */
+  if( instr[0]=='w' ? *wcstoimg_ptr : *imgtowcs_ptr )
+    error(EXIT_FAILURE, 0, "`%s' can only be called once.", optname);
+
+  /* Read the first token. */
+  c1=strtok_r(instr+9, delimiter, &saveptr);
+  if(c1==NULL)
+    error(EXIT_FAILURE, 0, "`%s' must end with `)'", optname);
+  if(*c1==')')
+    error(EXIT_FAILURE, 0, "`%s' has no defining column. Please specify "
+          "two column names or numbers with a slash between them for "
+          "example `%s(RA/DEC)'", optname, optname);
+
+  /* Read the second token and make sure it ends with `)'. Then set the `)'
+     as the string-NULL character. */
+  c2=strtok_r(NULL,        delimiter, &saveptr);
+  if(c2==NULL || c2[strlen(c2)-1]!=')')
+    error(EXIT_FAILURE, 0, "`%s' needs two input columns. Please "
+          "specify two columns with a slash between them for "
+          "example `%s(RA/DEC)'", optname, optname);
+  c2[strlen(c2)-1]='\0';
+
+  /* Make sure there aren't any more elements. */
+  c3=strtok_r(NULL, delimiter, &saveptr);
+  if(c3!=NULL)
+    error(EXIT_FAILURE, 0, "only two values can be given to `%s'", optname);
+
+  /* Add the column name/number to the list of columns to read, then keep
+     the pointer to the node (so we can later identify it). */
+  gal_list_str_add(new, c1, 1);
+  gal_list_str_add(new, c2, 1);
+
+  /* Keep this node's pointer. */
+  if(instr[0]=='w') *wcstoimg_ptr=*new;
+  else              *imgtowcs_ptr=*new;
+
+  /* Clean up (the `c1' and `c2') are actually within the allocated space
+     of `instr'. That is why we are setting the last argument of
+     `gal_list_str_add' to `1' (to allocate space for them). */
+  free(instr);
+}
+
+
+
+
+
 /* The columns can be given as comma-separated values to one option or
    multiple calls to the column option. Here, we'll check if the input list
    has comma-separated values. If they do then the list will be updated to
    be fully separate. */
 static void
-ui_columns_prepare(struct tableparams *p)
+ui_columns_prepare(struct tableparams *p, size_t *wcstoimg, size_t *imgtowcs)
 {
   size_t i;
   char **strarr;
   gal_data_t *strs;
   gal_list_str_t *tmp, *new=NULL;
+  gal_list_str_t *wcstoimg_ptr=NULL, *imgtowcs_ptr=NULL;
 
   /* Go over the whole original list (where each node may have more than
      one value separated by a comma. */
@@ -453,10 +532,18 @@ ui_columns_prepare(struct tableparams *p)
       strs=gal_options_parse_csv_strings_raw(tmp->v, NULL, 0);
       strarr=strs->array;
 
-      /* Go over all the elements and add them to the `new' list. */
+      /* Go over all the given colum names/numbers. */
       for(i=0;i<strs->size;++i)
         {
-          gal_list_str_add(&new, strarr[i], 0);
+          if(!strncmp(strarr[i],"wcstoimg(",9)
+             || !strncmp(strarr[i],"imgtowcs(",9))
+            ui_wcs_conversion(p, strarr[i], &wcstoimg_ptr, &imgtowcs_ptr, 
&new);
+          else
+            gal_list_str_add(&new, strarr[i], 0);
+
+          /* The pointer allocated string is either being used (and later
+             freed) else, or has already been freed. So its necessary to
+             set it to NULL. */
           strarr[i]=NULL;
         }
 
@@ -470,6 +557,16 @@ ui_columns_prepare(struct tableparams *p)
   /* Reverse the new list, then put it into `p->columns'. */
   gal_list_str_reverse(&new);
   p->columns=new;
+
+  /* If conversion is necessary, set the final column number. */
+  i=0;
+  if(wcstoimg_ptr || imgtowcs_ptr)
+    for(tmp=p->columns;tmp!=NULL;tmp=tmp->next)
+      {
+        if(wcstoimg_ptr==tmp) *wcstoimg=i-1;
+        if(imgtowcs_ptr==tmp) *imgtowcs=i-1;
+        ++i;
+      }
 }
 
 
@@ -769,9 +866,11 @@ static void
 ui_preparations(struct tableparams *p)
 {
   gal_list_str_t *lines;
+  size_t i, ncolnames, *colmatch;
   size_t nrange=0, origoutncols=0;
   struct gal_options_common_params *cp=&p->cp;
   size_t sortindout=GAL_BLANK_SIZE_T, *rangeindout=NULL;
+  size_t wcstoimg=GAL_BLANK_SIZE_T, imgtowcs=GAL_BLANK_SIZE_T;
 
   /* If there were no columns specified or the user has asked for
      information on the columns, we want the full set of columns. */
@@ -780,7 +879,7 @@ ui_preparations(struct tableparams *p)
 
 
   /* Prepare the column names. */
-  ui_columns_prepare(p);
+  ui_columns_prepare(p, &wcstoimg, &imgtowcs);
 
 
   /* If the input is from stdin, save it as `lines'. */
@@ -793,14 +892,30 @@ ui_preparations(struct tableparams *p)
                                &rangeindout);
 
 
+  /* If any conversions must be done, we need to know how many matches
+     there were for each column. */
+  ncolnames=gal_list_str_number(p->columns);
+  colmatch = ( (wcstoimg!=GAL_BLANK_SIZE_T || imgtowcs!=GAL_BLANK_SIZE_T)
+               ? gal_pointer_allocate(GAL_TYPE_SIZE_T, ncolnames, 1,
+                                      __func__, "colmatch")
+               : NULL );
+
+
   /* Read the necessary columns. */
   p->table=gal_table_read(p->filename, cp->hdu, lines, p->columns,
                           cp->searchin, cp->ignorecase, cp->minmapsize,
-                          NULL);
+                          colmatch);
   if(p->filename==NULL) p->filename="stdin";
   gal_list_str_free(lines, 1);
 
 
+  /* If we have a unit conversion, find the proper column to use. */
+  if(wcstoimg!=GAL_BLANK_SIZE_T)
+    { p->wcstoimg=0; for(i=0;i<wcstoimg;++i) p->wcstoimg += colmatch[i]; }
+  if(imgtowcs!=GAL_BLANK_SIZE_T)
+    { p->imgtowcs=0; for(i=0;i<imgtowcs;++i) p->imgtowcs += colmatch[i]; }
+
+
   /* If the range and sort options are requested, keep them as separate
      datasets. */
   if(p->range || p->sort)
@@ -838,6 +953,7 @@ ui_preparations(struct tableparams *p)
 
 
   /* Clean up. */
+  free(colmatch);
   if(rangeindout) free(rangeindout);
 }
 
diff --git a/bin/table/ui.h b/bin/table/ui.h
index 6247fe7..37f61a3 100644
--- a/bin/table/ui.h
+++ b/bin/table/ui.h
@@ -32,12 +32,14 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 /* Available letters for short options:
 
-   a b d e f g j k l m n p t u v w x y z
-   A B C E G H J L O Q R W X Y
+   a b d e f g j k l m n p t u v x y z
+   A B C E G H J L O Q R X Y
 */
 enum option_keys_enum
 {
   /* With short-option version. */
+  UI_KEY_WCSFILE         = 'w',
+  UI_KEY_WCSHDU          = 'W',
   UI_KEY_COLUMN          = 'c',
   UI_KEY_INFORMATION     = 'i',
   UI_KEY_COLINFOINSTDOUT = 'O',
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 0b797d4..13a0f5e 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -11729,6 +11729,47 @@ This option is not mandatory, if no specific columns 
are requested, all the
 input table columns are output. When this option is called multiple times,
 it is possible to output one column more than once.
 
+This option can also modify the input values based on certain
+operations. The operators are called just like a column name, allowing you
+to place the modified columns anywhere in the output table. For example the
+output of the two commands below (which will produce identical outputs)
+will be 5 columns. The first three columns are the input table's IDs and RA
+and Dec of each object. The fourth and fifth will be the pixel positions in
address@hidden that correspond to those positions.
+
address@hidden
+$ asttable table.fits -cID,RA,DEC,wcstoimg(RA/DEC) \
+           --wcsfile=image.fits
+$ asttable table.fits -cID,RA -cDEC -cwcstoimg(RA/DEC) \
+           --wcsfile=image.fits
address@hidden example
+
address@hidden @code
address@hidden wcstoimg(C1/C2)
+Assume that the column specfied by @code{C1} (either through name or
+number) is the RA and @code{C2} is the Declination. As output, the two
+outputs of this operator are the image coordinates. Note that within the
+operator, the delimiter/separator between column identifiers is a slash
+(address@hidden/}'). The WCS structure necessary for this conversion is derived
+from the @option{--wcshdu} extension/HDU of @option{--wcsfile}.
+
address@hidden imgtowcs(C1/C2)
+Similar to @code{wcstoimg}, except that the two input columns are assumed
+to be image coordinates and the output columns will be
address@hidden table
+
address@hidden -w STR
address@hidden --wcsfile=STR
+FITS file that contains the WCS to be used in the @code{wcstoimg} and
address@hidden operators of @option{--column} (see above). The extension
+name/number within the FITS file can be specified with @option{--wcshdu}.
+
address@hidden -W STR
address@hidden --wcshdu=STR
+FITS extension/HDU that contains the WCS to be used in the @code{wcstoimg}
+and @code{imgtowcs} operators of @option{--column} (see above). The FITS
+file name can be specified with @option{--wcsfile}.
+
 @item -O
 @itemx --colinfoinstdout
 @cindex Standard output