Fix LOGICAL type suport in TFITSIO

Solves a bug in the TFITSIO code when reading a FITS binary table that contained a column of type LOGICAL. The code in TFISTIO:LoadHDU only had tests for column types of strings or numeric so it was failing with the error “bad binary table datatype”
evheniyt · Feb 5, 2018 · 70433a5 · 70433a5
1 parent 55dff8c
commit 70433a5
Show file tree

Hide file tree

Showing 4 changed files with 561 additions and 36 deletions.
diff --git a/.gitattributes b/.gitattributes
@@ -17,3 +17,4 @@
 *.jpg binary
 *.pdf binary
 *.root binary
+*.fits binary
diff --git a/graf2d/fitsio/src/TFITS.cxx b/graf2d/fitsio/src/TFITS.cxx
@@ -85,6 +85,7 @@ several methods to manage them.
 #include "fitsio.h"
 #include <stdlib.h>
 
+
 ClassImp(TFITSHDU);
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -380,6 +381,7 @@ Bool_t TFITSHDU::LoadHDU(TString& filepath_filter)
 
       for (colnum = 0, cellindex = 0; colnum < fNColumns; colnum++) {
          fits_get_coltype(fp, colnum+1, &typecode, &repeat, &width, &status);
+
          if (status) goto ERR;
 
          if ((typecode == TDOUBLE) || (typecode == TSHORT) || (typecode == TLONG)
@@ -441,60 +443,89 @@ Bool_t TFITSHDU::LoadHDU(TString& filepath_filter)
 
                fColumnsInfo[colnum].fDim = (Int_t) repeat;
 
-               double *array;
+               double *array = 0;
+               char *arrayl = 0;
 
                if (repeat > 0) {
-                  array = new double [table_rows * repeat]; //Hope you got a big machine! Ask China otherwise :-)
-                  fits_read_col(fp, TDOUBLE, colnum+1, 1, 1, table_rows * repeat, &nulval, array, &anynul, &status);
 
-                  if (status) {
-                     delete [] array;
-                     goto ERR;
+                  if (typecode == TLOGICAL) {
+                     arrayl = new char[table_rows * repeat];
+                     fits_read_col(fp, TLOGICAL, colnum + 1, 1, 1, table_rows * repeat, &nulval, arrayl, &anynul,
+                                   &status);
+                     if (status) {
+                        delete[] arrayl;
+                        goto ERR;
+                     }
+                  } else {
+                     array = new double[table_rows * repeat]; // Hope you got a big machine! Ask China otherwise :-)
+                     fits_read_col(fp, TDOUBLE, colnum + 1, 1, 1, table_rows * repeat, &nulval, array, &anynul,
+                                   &status);
+                     if (status) {
+                        delete[] array;
+                        goto ERR;
+                     }
                   }
+
                } else {
-                  //No elements: set dummy
-                  array = new double [table_rows];
+                  // No elements: set dummy
+                  array = new double[table_rows];
                   for (long row = 0; row < table_rows; row++) {
                      array[row] = 0.0;
                   }
                }
 
-               //Save values
+               // Save values
                if (repeat == 1) {
-                  //Scalar
-                  for (long row = 0; row < table_rows; row++) {
-                     fCells[cellindex++].fRealNumber = array[row];
+                  // Scalar
+                  if (typecode == TLOGICAL) {
+                     for (long row = 0; row < table_rows; row++) {
+                        int temp = (signed char)arrayl[row];
+                        fCells[cellindex++].fRealNumber = (double)temp;
+                     }
+                     delete[] arrayl;
+                  } else {
+                     for (long row = 0; row < table_rows; row++) {
+                        fCells[cellindex++].fRealNumber = array[row];
+                     }
+                     delete[] array;
                   }
                } else if (repeat > 1) {
-                  //Vector
-                  for (long row = 0; row < table_rows; row++) {
-                     double *vec = new double [repeat];
-                     long offset = row * repeat;
-                     for (long component = 0; component < repeat; component++) {
-                        vec[component] = array[offset++];
+                  // Vector
+                  if (typecode == TLOGICAL) {
+                     for (long row = 0; row < table_rows; row++) {
+                        double *vec = new double[repeat];
+                        long offset = row * repeat;
+                        for (long component = 0; component < repeat; component++) {
+                           int temp = (signed char)arrayl[offset++];
+                           vec[component] = (double)temp;
+                        }
+                        fCells[cellindex++].fRealVector = vec;
+                     }
+                     delete[] arrayl;
+                  } else {
+                     for (long row = 0; row < table_rows; row++) {
+                        double *vec = new double[repeat];
+                        long offset = row * repeat;
+                        for (long component = 0; component < repeat; component++) {
+                           vec[component] = array[offset++];
+                        }
+                        fCells[cellindex++].fRealVector = vec;
                      }
-                     fCells[cellindex++].fRealVector = vec;
+                     delete[] array;
                   }
                }
-
-               delete [] array;
-
             }
-
          } else {
             Warning("LoadHDU", "error opening FITS file. Column type %d is currently not supported", typecode);
          }
       }
 
-
-
       if (hdutype == ASCII_TBL) {
-         //ASCII table
+         // ASCII table
 
       } else {
-         //Binary table
+         // Binary table
       }
-
    }
 
    // Close file
@@ -1058,13 +1089,16 @@ TVectorD* TFITSHDU::GetArrayRow(UInt_t row)
    }
 
    offset = W * row;
-   TVectorD *vec = new TVectorD(W);
-   double *v = vec->GetMatrixArray();
+   double *v = new double[W];
 
    for (i = 0; i < W; i++) {
       v[i] = fPixels->GetAt(offset+i);
    }
 
+   TVectorD *vec = new TVectorD(W, v);
+
+   delete [] v;
+
    return vec;
 }
 
@@ -1094,13 +1128,15 @@ TVectorD* TFITSHDU::GetArrayColumn(UInt_t col)
       return 0;
    }
 
-   TVectorD *vec = new TVectorD(H);
-   double *v = vec->GetMatrixArray();
+   double *v = new double[H];
 
    for (i = 0; i < H; i++) {
       v[i] = fPixels->GetAt(W*i+col);
    }
 
+   TVectorD *vec = new TVectorD(H, v);
+
+   delete [] v;
 
    return vec;
 }
@@ -1208,13 +1244,15 @@ TVectorD* TFITSHDU::GetTabRealVectorColumn(Int_t colnum)
 
    Int_t offset = colnum * fNRows;
 
-   TVectorD *res = new TVectorD(fNRows);
-   Double_t *arr = res->GetMatrixArray();
+   Double_t *arr = new Double_t[fNRows];
 
    for (Int_t row = 0; row < fNRows; row++) {
       arr[row] = fCells[offset + row].fRealNumber;
    }
 
+   TVectorD *res = new TVectorD();
+   res->Use(fNRows, arr);
+
    return res;
 }
 
@@ -1245,13 +1283,15 @@ TVectorD* TFITSHDU::GetTabRealVectorColumn(const char *colname)
 
    Int_t offset = colnum * fNRows;
 
-   TVectorD *res = new TVectorD(fNRows);
-   Double_t *arr = res->GetMatrixArray();
+   Double_t *arr = new Double_t[fNRows];
 
    for (Int_t row = 0; row < fNRows; row++) {
       arr[row] = fCells[offset + row].fRealNumber;
    }
 
+   TVectorD *res = new TVectorD();
+   res->Use(fNRows, arr);
+
    return res;
 }
 

diff --git a/tutorials/fitsio/FITS_tutorial7.C b/tutorials/fitsio/FITS_tutorial7.C
@@ -0,0 +1,63 @@
+/// \file
+/// \ingroup tutorial_FITS
+/// Open a FITS file that contains a catalog of astronomical objects
+/// and dump some of its columns
+///
+/// \macro_code
+///
+/// \author Elizabeth Buckley-Geer
+
+void FITS_tutorial7()
+{
+
+   printf("\n\n--------------------------------\n");
+   printf("WELCOME TO FITS tutorial #7 !!!!\n");
+   printf("--------------------------------\n");
+   printf("We are going to open a table from a FITS file\n");
+   printf("and print out three columns for some of the objects.\n");
+   printf("This table contains a logical data type so this tutorial tests\n");
+   printf("that we can read it correctly\n\n");
+
+   TString dir = gSystem->DirName(__FILE__);
+
+   // Open the table
+   TFITSHDU *hdu = new TFITSHDU(dir + "/sample5.fits[1]");
+   if (hdu == 0) {
+      printf("ERROR: could not access the HDU\n");
+      return;
+   }
+
+   TVectorD *vec1;
+   TVectorD *vec2;
+   TVectorD *vec3;
+   TVectorD *vec4;
+
+   // Read the ra, dec, flux_g and brick_primary columns
+
+   vec1 = hdu->GetTabRealVectorColumn("ra");
+   vec2 = hdu->GetTabRealVectorColumn("dec");
+   vec3 = hdu->GetTabRealVectorColumn("flux_g");
+   vec4 = hdu->GetTabRealVectorColumn("brick_primary");
+
+   Double_t gflux, ra, dec, bp;
+
+   for (Int_t i = vec1->GetLwb(); i <= vec1->GetUpb(); i++) {
+
+      bp = (*vec4)[i];
+      gflux = (*vec3)[i];
+      ra = (*vec1)[i];
+      dec = (*vec2)[i];
+
+      if (bp) {
+         printf("RA %f DEC %f G-FLUX %f\n", ra, dec, gflux);
+      }
+   }
+
+   // Clean up
+
+   delete vec1;
+   delete vec2;
+   delete vec3;
+   delete vec4;
+   delete hdu;
+}
diff --git a/tutorials/fitsio/sample5.fits b/tutorials/fitsio/sample5.fits