[gnuastro-commits] master d09e850 04/16: FITS table column information printed by Table

[Mohammad Akhlaghi]

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

    FITS table column information printed by Table
    
    Table can now print the basic column information of a given FITS ASCII or
    Binary table (so far only tested on a binary table). This is far from
    complete for this utility, but for the time being it is a significant step
    (it actually prints something useful!).
    
    The type of the `gal_fits_read_hdu' function was changed from the exact
    CFITSIO macro to an image or table binary value so we could use this
    function in reading tables too.
---
 Makefile.am           |    2 +-
 lib/fits.c            |   99 ++++++++++++++++++++++++++++---
 lib/gnuastro/fits.h   |   12 +++-
 src/imgcrop/imgcrop.c |    4 +-
 src/imgcrop/ui.c      |    2 +-
 src/table/args.h      |   32 +++++-----
 src/table/main.h      |   19 +++++-
 src/table/table.c     |   76 +++++++++++++++++++++++-
 src/table/ui.c        |  157 ++++++++++++++++++++++++++++++++++++++++++++-----
 9 files changed, 353 insertions(+), 50 deletions(-)

diff --git a/Makefile.am b/Makefile.am
index fae1bf4..1c819fa 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -108,7 +108,7 @@ SUBDIRS = bootstrapped/lib $(MAYBE_GNULIBCHECK) lib 
$(MAYBE_ARITHMETIC)     \
 $(MAYBE_CONVERTT) $(MAYBE_CONVOLVE) $(MAYBE_COSMICCAL) $(MAYBE_HEADER)      \
 $(MAYBE_IMGCROP) $(MAYBE_IMGSTAT) $(MAYBE_IMGWARP) $(MAYBE_MKCATALOG)       \
 $(MAYBE_MKNOISE) $(MAYBE_MKPROF) $(MAYBE_NOISECHISEL) $(MAYBE_SUBTRACTSKY)  \
-$(MAYBE_TEMPLATE) doc tests
+$(MAYBE_TABLE) $(MAYBE_TEMPLATE) doc tests
 
 
 
diff --git a/lib/fits.c b/lib/fits.c
index 1168bd3..922ba9a 100644
--- a/lib/fits.c
+++ b/lib/fits.c
@@ -733,7 +733,7 @@ hdutypestring(int hdutype)
 /* Check the desired HDU in a FITS image and also if it has the
    desired type. */
 void
-gal_fits_read_hdu(char *filename, char *hdu, int desiredtype,
+gal_fits_read_hdu(char *filename, char *hdu, unsigned char img0_tab1,
                   fitsfile **outfptr)
 {
   size_t len;
@@ -758,11 +758,24 @@ gal_fits_read_hdu(char *filename, char *hdu, int 
desiredtype,
   if (fits_get_hdu_type(fptr, &hdutype, &status) )
     gal_fits_io_error(status, NULL);
 
-  if(hdutype!=desiredtype)
-    error(EXIT_FAILURE, 0, "%s: HDU %s is %s, not %s",
-          filename, hdu, hdutypestring(hdutype),
-          hdutypestring(desiredtype));
 
+  /* Check if the type of the HDU is the expected type. We could have
+     written these as && conditions, but this is easier to read, it makes
+     no meaningful difference to the compiler. */
+  if(img0_tab1)
+    {
+      if(hdutype==IMAGE_HDU)
+        error(EXIT_FAILURE, 0, "%s: HDU %s is an image, not a table",
+              filename, hdu);
+    }
+  else
+    {
+      if(hdutype!=IMAGE_HDU)
+        error(EXIT_FAILURE, 0, "%s: HDU %s is %s, not an image",
+              filename, hdu, hdutypestring(hdutype));
+    }
+
+  /* Clean up. */
   free(ffname);
 }
 
@@ -1309,7 +1322,7 @@ gal_fits_read_wcs(char *filename, char *hdu, size_t 
hstartwcs,
   fitsfile *fptr;
 
   /* Check HDU for realistic conditions: */
-  gal_fits_read_hdu(filename, hdu, IMAGE_HDU, &fptr);
+  gal_fits_read_hdu(filename, hdu, 0, &fptr);
 
   /* Read the WCS information: */
   gal_fits_read_wcs_from_pointer(fptr, nwcs, wcs, hstartwcs, hendwcs);
@@ -1332,7 +1345,7 @@ gal_fits_read_wcs(char *filename, char *hdu, size_t 
hstartwcs,
    the macros in fitsarrayvv.h and depends on the type of the data.*/
 int
 gal_fits_hdu_to_array(char *filename, char *hdu, int *bitpix,
-                                void **array, size_t *s0, size_t *s1)
+                      void **array, size_t *s0, size_t *s1)
 {
   void *bitblank;
   fitsfile *fptr;
@@ -1340,7 +1353,7 @@ gal_fits_hdu_to_array(char *filename, char *hdu, int 
*bitpix,
   long naxes[2], fpixel[]={1,1};
 
   /* Check HDU for realistic conditions: */
-  gal_fits_read_hdu(filename, hdu, IMAGE_HDU, &fptr);
+  gal_fits_read_hdu(filename, hdu, 0, &fptr);
 
   /* Get the bitpix and size of the image: */
   gal_fits_img_bitpix_size(fptr, bitpix, naxes);
@@ -1516,6 +1529,76 @@ gal_fits_atof_correct_wcs(char *filename, char *hdu, int 
bitpix,
 
 
 /**************************************************************/
+/**********                 Table                  ************/
+/**************************************************************/
+/* The values to the TFORM header keyword are single letter capital
+   letters, but that is useless in identifying the data type of the
+   column. So this function will do the conversion based on the CFITSIO
+   manual.*/
+int
+gal_fits_tform_to_dtype(char tform)
+{
+  switch(tform)
+    {
+    case 'X':
+      return TBIT;
+    case 'B':
+      return TBYTE;
+    case 'L':
+      return TLOGICAL;
+    case 'A':
+      return TSTRING;
+    case 'I':
+      return TSHORT;
+    case 'J':
+      return TLONG;
+    case 'K':
+      return TLONGLONG;
+    case 'E':
+      return TFLOAT;
+    case 'D':
+      return TDOUBLE;
+    case 'C':
+      return TCOMPLEX;
+    case 'M':
+      return TDBLCOMPLEX;
+    case 'S':
+      return TSBYTE;
+    case 'V':
+      return TUINT;
+    case 'U':
+      return TUSHORT;
+    default:
+      error(EXIT_FAILURE, 0, "'%c' is not a recognized CFITSIO value for "
+            "the TFORMn header keyword(s).", tform);
+    }
+
+  error(EXIT_FAILURE, 0, "A bug! Please contact us so we can fix this. "
+        "For some reason, control has reached to the end of the "
+        "gal_fits_tform_to_dtype function in fits.c.");
+  return -1;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/**************************************************************/
 /**********          Check prepare file            ************/
 /**************************************************************/
 /* We have the name of the input file. But in most cases, the files
diff --git a/lib/gnuastro/fits.h b/lib/gnuastro/fits.h
index a0456a4..216fdb4 100644
--- a/lib/gnuastro/fits.h
+++ b/lib/gnuastro/fits.h
@@ -172,7 +172,7 @@ void
 gal_fits_img_bitpix_size(fitsfile *fptr, int *bitpix, long *naxis);
 
 void
-gal_fits_read_hdu(char *filename, char *hdu, int desiredtype,
+gal_fits_read_hdu(char *filename, char *hdu, unsigned char img0_tab1,
                   fitsfile **outfptr);
 
 void *
@@ -213,6 +213,16 @@ gal_fits_atof_correct_wcs(char *filename, char *hdu, int 
bitpix,
 
 
 /**************************************************************/
+/**********                  Table                 ************/
+/**************************************************************/
+int
+gal_fits_tform_to_dtype(char tform);
+
+
+
+
+
+/**************************************************************/
 /**********          Check prepare file            ************/
 /**************************************************************/
 void
diff --git a/src/imgcrop/imgcrop.c b/src/imgcrop/imgcrop.c
index 849ed75..2f17b09 100644
--- a/src/imgcrop/imgcrop.c
+++ b/src/imgcrop/imgcrop.c
@@ -90,7 +90,7 @@ imgmodecrop(void *inparam)
   /* The whole catalog is from one image, so you can get the
      information here:*/
   img=&p->imgs[crp->imgindex];
-  gal_fits_read_hdu(img->name, cp->hdu, IMAGE_HDU, &crp->infits);
+  gal_fits_read_hdu(img->name, cp->hdu, 0, &crp->infits);
 
   /* Go over all the outputs that are assigned to this thread: */
   for(i=0;crp->indexs[i]!=GAL_THREADS_NON_THRD_INDEX;++i)
@@ -181,7 +181,7 @@ wcsmodecrop(void *inparam)
         if(radecoverlap(crp))
           {
             gal_fits_read_hdu(p->imgs[crp->imgindex].name, p->cp.hdu,
-                                        IMAGE_HDU, &crp->infits);
+                              0, &crp->infits);
 
             if(log->name==NULL) cropname(crp);
 
diff --git a/src/imgcrop/ui.c b/src/imgcrop/ui.c
index 82679ff..c10e2ac 100644
--- a/src/imgcrop/ui.c
+++ b/src/imgcrop/ui.c
@@ -621,7 +621,7 @@ preparearrays(struct imgcropparams *p)
       status=0;
       img=&p->imgs[i];
       gal_linkedlist_pop_from_stll(&p->up.gal_linkedlist_stll, &img->name);
-      gal_fits_read_hdu(img->name, p->cp.hdu, IMAGE_HDU, &tmpfits);
+      gal_fits_read_hdu(img->name, p->cp.hdu, 0, &tmpfits);
       gal_fits_img_bitpix_size(tmpfits, &p->bitpix, img->naxes);
       gal_fits_read_wcs_from_pointer(tmpfits, &img->nwcs, &img->wcs,
                                      p->hstartwcs, p->hendwcs);
diff --git a/src/table/args.h b/src/table/args.h
index 74fa061..0583dcd 100644
--- a/src/table/args.h
+++ b/src/table/args.h
@@ -155,39 +155,35 @@ parse_opt(int key, char *arg, struct argp_state *state)
     /* Operating modes: */
     case 'i':
       p->information=1;
-      p->informationset=1;
+      p->up.informationset=1;
       break;
 
 
     /* Read the non-option arguments: */
     case ARGP_KEY_ARG:
 
-      /* See what type of input value it is and put it in. */
-      if( gal_fits_name_is_fits(arg) )
-        {
-          if(p->up.inputname)
-            argp_error(state, "only one input image should be given");
-          else
-            p->up.inputname=arg;
-        }
-      else
-        argp_error(state, "%s is not a valid file type", arg);
-      break;
+      /* Table gets only one input argument. */
+      if(p->up.inputset)
+        argp_error(state, "only one input file should be given");
 
+      /* This is the first (and must be only) argument. */
+      p->up.inputset=1;
 
+      /* See what type of input this is, and save the value. */
+      if( gal_fits_name_is_fits(arg) )
+        p->up.fitsname=arg;
+      else
+        p->up.txtname=arg;
 
+      break;
 
 
     /* The command line options and arguments are finished. */
     case ARGP_KEY_END:
       if(p->cp.setdirconf==0 && p->cp.setusrconf==0
          && p->cp.printparams==0)
-        {
-          if(state->arg_num==0)
-            argp_error(state, "no argument given");
-          if(p->up.inputname==NULL)
-            argp_error(state, "no input FITS image(s) provided");
-        }
+        if(state->arg_num==0)
+          argp_error(state, "no argument given");
       break;
 
 
diff --git a/src/table/main.h b/src/table/main.h
index 0f3c1fb..7a15c05 100644
--- a/src/table/main.h
+++ b/src/table/main.h
@@ -36,9 +36,16 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 
 
+
+
+
+/* User interface structure. */
 struct uiparams
 {
-  char             *inputname;  /* Name of input file.             */
+  char              *fitsname;  /* Name of input FITS file.        */
+  char               *txtname;  /* Name of input text file.        */
+
+  int                inputset;
   int          informationset;
 };
 
@@ -46,6 +53,7 @@ struct uiparams
 
 
 
+/* Main program parameters structure */
 struct tableparams
 {
   /* Other structures: */
@@ -54,6 +62,15 @@ struct tableparams
 
   /* Input: */
   int             information;  /* ==1: only print information.       */
+  fitsfile           *fitsptr;  /* FITS pointer (input or output).    */
+
+  /* FITS table */
+  size_t                nrows;  /* Number of rows in table.           */
+  size_t                ncols;  /* Number of columns in table.        */
+  int               *typecode;  /* Type of data in column.            */
+  char                **tform;  /* TFORM (another format for type).   */
+  char                **ttype;  /* Column name (one word).            */
+  char                **tunit;  /* Unit of values in column.          */
 
   /* Internal: */
   int                onlyview;
diff --git a/src/table/table.c b/src/table/table.c
index ab560f2..eb742f5 100644
--- a/src/table/table.c
+++ b/src/table/table.c
@@ -29,14 +29,86 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <stdlib.h>
 
 #include <gnuastro/fits.h>
-#include <gnuastro/linkedlist.h>
 
 #include "main.h"
 
 
 
+
+/* Print all the column information: */
 void
-table(struct tableparams *p)
+printinfo(struct tableparams *p)
 {
+  size_t i;
+  char *typestring=NULL;
+
+  printf("%s (hdu: %s)\n", p->up.fitsname, p->cp.hdu);
+  printf("Number of rows: %lu\n", p->nrows);
+  printf("Column information\n");
+  printf("------------------\n");
+  for(i=1;i<=p->ncols;++i)
+    {
+      switch(p->typecode[i])
+        {
+        case TBIT:
+          typestring="bit";
+          break;
+        case TBYTE:
+          typestring="byte";
+          break;
+        case TLOGICAL:
+          typestring="logicals";
+          break;
+        case TSTRING:
+          typestring="string";
+          break;
+        case TSHORT:
+          typestring="short";
+          break;
+        case TLONG:
+          typestring="long";
+          break;
+        case TLONGLONG:
+          typestring="longlong";
+          break;
+        case TFLOAT:
+          typestring="float";
+          break;
+        case TDOUBLE:
+          typestring="double";
+          break;
+        case TCOMPLEX:
+          typestring="complex";
+          break;
+        case TDBLCOMPLEX:
+          typestring="dblcomplex";
+          break;
+        case TSBYTE:
+          typestring="signed byte";
+          break;
+        case TUINT:
+          typestring="unsigned int";
+          break;
+        case TUSHORT:
+          typestring="unsigned short";
+          break;
+        default:
+          error(EXIT_FAILURE, 0, "%d (from TFORM%lu='%c') is not a "
+                "recognized CFITSIO datatype.",
+                p->typecode[i], i, p->tform[i][0]);
+        }
+      printf("%-3lu %-25s %-20s %s\n", i, p->ttype[i], typestring,
+             p->tunit[i]);
+    }
+
+}
+
 
+
+
+/* Top level function */
+void
+table(struct tableparams *p)
+{
+  if(p->information) printinfo(p);
 }
diff --git a/src/table/ui.c b/src/table/ui.c
index 4903278..bdd75f4 100644
--- a/src/table/ui.c
+++ b/src/table/ui.c
@@ -102,30 +102,23 @@ readconfig(char *filename, struct tableparams *p)
 
       /* Inputs: */
       if(strcmp(name, "hdu")==0)
-        {
-          if(cp->hduset) continue;
-          errno=0;
-          cp->hdu=malloc(strlen(value)+1);
-          if(cp->hdu==NULL)
-            error(EXIT_FAILURE, errno, "space for HDU");
-          strcpy(cp->hdu, value);
-          cp->hduset=1;
-        }
+        gal_checkset_allocate_copy_set(value, &cp->hdu, &cp->outputset);
 
 
 
       /* Outputs */
-
+      else if(strcmp(name, "output")==0)
+        gal_checkset_allocate_copy_set(value, &cp->output, &cp->outputset);
 
 
 
       /* Operating modes: */
       else if (strcmp(name, "information")==0)
         {
-          if(up->informationset) continue;
+          if(p->up.informationset) continue;
           gal_checkset_int_zero_or_one(value, &p->information, name,
                                        key, SPACK, filename, lineno);
-          up->informationset=1;
+          p->up.informationset=1;
         }
 
 
@@ -210,10 +203,138 @@ checkifset(struct tableparams *p)
 /**************************************************************/
 /***************       Sanity Check         *******************/
 /**************************************************************/
+/* This function will read all the table information from a FITS table HDU
+   and store them in arrays for use later. It is mainly good for getting
+   all the information in a FITS table HDU. This function will only go
+   through the header keywords once and does not depend on the ordering of
+   the keywords, so it is much more efficient than having to ask for each
+   column's information separately.*/
+void
+readallcolinfo(char *filename, char *hdu, size_t *nrows,
+               size_t *ncols, int **otypecode, char ***otform,
+               char ***ottype, char ***otunit, fitsfile **outfptr)
+{
+  long lnrows;
+  size_t index;
+  char *tailptr, **c, **fc;
+  char **tform, **ttype, **tunit;
+  int i, status=0, incols, *typecode;
+  char keyname[FLEN_KEYWORD]="XXXXXXXXXXXXX", value[FLEN_VALUE];
+
+  /* Set the FITS pointer and check the HDU type. */
+  gal_fits_read_hdu(filename, hdu, 1, outfptr);
+
+  /* Get the size of the table. incols is incremented because the indexs of
+     the FITS keywords start from 1, not zero. */
+  fits_get_num_rows(*outfptr, &lnrows, &status);
+  fits_get_num_cols(*outfptr, &incols, &status);
+  *ncols=incols++;
+  *nrows=lnrows;
+
+  /* Allocate the arrays to keep the column information. Initialize the
+     arrays with a NULL pointer to make sure that they are all found in the
+     end if they are necessary.*/
+  errno=0; typecode=*otypecode=malloc(incols * sizeof *typecode);
+  if(typecode==NULL)
+    error(EXIT_FAILURE, errno, "%lu bytes for typecode",
+          incols * sizeof *typecode);
+  errno=0; tform=*otform=malloc(incols * sizeof *tform);
+  if(tform==NULL)
+    error(EXIT_FAILURE, errno, "%lu bytes for tform",
+          incols * sizeof *tform);
+  fc=(c=tform)+incols; do *c++=NULL; while(c<fc);
+  errno=0; ttype=*ottype=malloc(incols * sizeof *ttype);
+  if(ttype==NULL)
+    error(EXIT_FAILURE, errno, "%lu bytes for ttype",
+          incols * sizeof *ttype);
+  fc=(c=ttype)+incols; do *c++=NULL; while(c<fc);
+  errno=0; tunit=*otunit=malloc(incols * sizeof *tunit);
+  if(tunit==NULL)
+    error(EXIT_FAILURE, errno, "%lu bytes for tunit",
+          incols * sizeof *tunit);
+  fc=(c=tunit)+incols; do *c++=NULL; while(c<fc);
+
+  /* Read all the keywords one by one and if they match, then put them in
+     the correct value. Note that we are starting from keyword 9 because
+     according to the FITS standard, the first 8 keys in a FITS table are
+     reserved. */
+  for(i=9; strcmp(keyname, "END"); ++i)
+    {
+      /* Read the next keyword. */
+      fits_read_keyn(*outfptr, i, keyname, value, NULL, &status);
+
+      /* Check the type of the keyword. */
+      if(strncmp(keyname, "TFORM", 5)==0)
+        {
+          /* The values to TFORM are only a single character, so start the
+             pointer to copy at 1 and put the string terminator at 3. */
+          value[2]='\0';
+          index=strtoul(&keyname[5], &tailptr, 10);
+          if(index<incols)
+            {
+              gal_checkset_allocate_copy(&value[1], &tform[index] );
+              typecode[index]=gal_fits_tform_to_dtype(value[1]);
+            }
+        }
+      else if(strncmp(keyname, "TTYPE", 5)==0)
+        {
+          /* All strings in CFITSIO start and finish with single quotation
+             marks, CFITSIO puts them in itsself, so if we don't remove
+             them here, we might have duplicates later, its easier to just
+             remove them to have a simple string that might be used else
+             where too (without the single quotes). */
+          value[strlen(value)-1]='\0';
+          index=strtoul(&keyname[5], &tailptr, 10);
+          if(index<incols)
+            gal_checkset_allocate_copy(&value[1], &ttype[index] );
+        }
+      else if(strncmp(keyname, "TUNIT", 5)==0)
+        {
+          /* similar to ttype, see above.*/
+          value[strlen(value)-1]='\0';
+          index=strtoul(&keyname[5], &tailptr, 10);
+          if(index<incols)
+            gal_checkset_allocate_copy(&value[1], &tunit[index] );
+        }
+    }
+
+  /* Check if all values were set: */
+  for(i=1;i<=*ncols;++i)
+    {
+      if(!tform[i])
+         error(EXIT_FAILURE, 0, "TFORM%d could not be found in header", i);
+      if(!ttype[i])
+         error(EXIT_FAILURE, 0, "TTYPE%d could not be found in header", i);
+      if(!tunit[i])
+         error(EXIT_FAILURE, 0, "TUNIT%d could not be found in header", i);
+    }
+
+
+  /* For a checkup:
+  for(i=1;i<=*ncols;++i)
+    printf("%d: %s, %s, %s\n", i, tform[i], ttype[i], tunit[i]);
+  */
+
+}
+
+
+
+
+
 void
 sanitycheck(struct tableparams *p)
 {
 
+  /* Set the FITS pointer and check the type of the fits file. */
+  if(p->up.fitsname)
+    readallcolinfo(p->up.fitsname, p->cp.hdu, &p->nrows,
+                   &p->ncols, &p->typecode, &p->tform,
+                   &p->ttype, &p->tunit, &p->fitsptr);
+  else
+    error(EXIT_FAILURE, 0, "Table is a new addition to Gnuastro and "
+          "and under heavy development, it currently doesn't support "
+          "anything other than a FITS binary table.");
+
 }
 
 
@@ -267,6 +388,7 @@ preparearrays(struct tableparams *p)
 void
 setparams(int argc, char *argv[], struct tableparams *p)
 {
+  struct uiparams *up=&p->up;
   struct gal_commonparams *cp=&p->cp;
 
   /* Set the non-zero initial values, the structure was initialized to
@@ -276,6 +398,9 @@ setparams(int argc, char *argv[], struct tableparams *p)
   cp->numthreads    = num_processors(NPROC_CURRENT);
   cp->removedirinfo = 1;
 
+  /* Initialize this utility's special variables. */
+  up->txtname=up->fitsname=NULL;
+
   /* Read the arguments. */
   errno=0;
   if(argp_parse(&thisargp, argc, argv, 0, 0, p))
@@ -323,13 +448,13 @@ setparams(int argc, char *argv[], struct tableparams *p)
 void
 freeandreport(struct tableparams *p)
 {
+  int status=0;
+
   /* Free the allocated arrays: */
   free(p->cp.hdu);
   free(p->cp.output);
 
-  /* Close the FITS file:
-  if(fits_close_file(p->fptr, &status))
+  /* Close the FITS file: */
+  if(p->up.fitsname && fits_close_file(p->fitsptr, &status))
     gal_fits_io_error(status, NULL);
-  */
-
 }