Changeset 32228

Timestamp:

Aug 30, 2011, 3:47:39 PM (15 years ago)

Author:

bills

Message:

working version of psFitsReadTableNew and psFitsWriteTableNew

Location:

trunk/psLib/src/fits

Files:

• Makefile.am (modified) (2 diffs)
• psFitsTableNew.c (modified) (14 diffs)
• psFitsTableNew.h (modified) (2 diffs)

trunk/psLib/src/fits/Makefile.am

@@ -9 +9 @@
         psFitsImage.c \
         psFitsTable.c \
+        psFitsTableNew.c \
         psFitsFloat.c \
         psFitsFloatFile.c \
@@ -20 +21 @@
         psFitsImage.h \
         psFitsTable.h \
+        psFitsTableNew.h \
         psFitsFloat.h \
         psFitsFloatFile.h \

trunk/psLib/src/fits/psFitsTableNew.c

@@ -96 +96 @@
 }
 
+void breakhere()
+{
+}
+
+void
+freeTable(psFitsTable *table) {
+    for (int col = 0; col < table->numCols; col++) {
+        psFitsTableColumn *column = &table->columns[col];
+        psFree(column->name);
+        // all of the members in the data union are pointers so just pick S32
+        psFree(column->data.S32);
+    }
+    psFree(table->columns);
+    psFree(table->index);
+}
+
+psFitsTable* psFitsTableAlloc(int numCols, int numRows) {
+    if (numCols <= 0) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                "Can't allocate a psFitsTable with %d columns.", numCols);
+        return NULL;
+    }
+
+    psFitsTable *table = psAlloc(sizeof(psFitsTable));
+    table->index = psMetadataAlloc();
+    table->numRows = numRows;
+    table->numCols = numCols;
+    table->columns = psAlloc(sizeof(psFitsTableColumn) * numCols);
+    memset(table->columns, sizeof(psFitsTableColumn) * numCols, 0);
+    psMemSetDeallocator(table, (psFreeFunc) freeTable);
+
+    return table;
+}
 
 psFitsTable* psFitsReadTableNew(const psFits* fits)
@@ -122 +155 @@
 
     psTrace("psLib.fits",5,"Table size is %ix%li\n",numCols, numRows);
-    // the row parameter in the proper range?
-#ifdef notdef
-    if (row < 0 || row >= numRows) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
-                _("Specified row, %d, is not valid for current table of %ld rows."),
-                row, numRows);
-        return NULL;
-    }
-#endif
-
-    psFitsTable *table = psAlloc(sizeof(psFitsTable));
-    table->index = psMetadataAlloc(); \
-    table->numRows = numRows;
-    table->numCols = numCols;
-    table->columns = psAlloc(sizeof(psFitsTableColumn) * numCols);
-//    psMetadata* data = psMetadataAlloc();
+
+    psFitsTable *table = psFitsTableAlloc(numCols, numRows);
 
     int hdutype;                        // Type of HDU: need to distinguish ASCII and binary tables
@@ -151 +170 @@
     for (int col = 1; col <= numCols; col++) {
         // get the column name
+        double tscal;
+        double tzero;
         if (hdutype == BINARY_TBL) {
             fits_get_bcolparms(fits->fd, col, name,
-                               NULL, NULL, NULL, NULL, NULL, NULL, NULL, &status);
+                               NULL, NULL, NULL, &tscal, &tzero, NULL, NULL, &status);
         } else {
             fits_get_acolparms(fits->fd, col, name,
@@ -170 +191 @@
                 column->data.TYPE = psAlloc(sizeof(ps##TYPE)*table->numRows); \
                 if (repeat == 1) { \
+                    column->elementSize = 1; \
                     NATIVETYPE *values = (NATIVETYPE *) psAlloc(sizeof(NATIVETYPE) * table->numRows); \
                     int anynul = 0; \
+                    breakhere(); \
                     fits_read_col(fits->fd, FITSTYPE, col, 1, \
                                   1, numRows, NULL, values, &anynul, &status); \
@@ -180 +203 @@
                             col, name, typecode, repeat); \
                 } else { \
+                    column->elementSize = repeat; \
                     column->type = PS_DATA_VECTOR; \
+                    column->vectorType = PS_DATA_##TYPE; \
                     column->data.vec = psAlloc(sizeof(psVector*) * table->numRows); \
                     for (int row = 0; row < table->numRows; row++) { \
@@ -192 +217 @@
                         } \
                         column->data.vec[row] = vec; \
-                        /* psMetadataAdd(data,PS_LIST_TAIL, name, PS_DATA_VECTOR, "", vec); */ \
                         psFree(values); \
-                        /* psFree(vec); */ \
                     } \
                     \
@@ -201 +224 @@
             }
 
+        // Handle cases where the data is "really" unsigned
+        if (typecode == TLONG && tzero == 2147483648) {
+            typecode = TULONG;
+        } else if (typecode == TSHORT && tzero == 32768) {
+            typecode = TUSHORT;
+        }
+
+        // Save the column number for this column in the index metadata
         psMetadataAddS32(table->index, PS_LIST_TAIL, name, 0, NULL, col-1);
+
         psFitsTableColumn *column = &table->columns[col-1];
-        // column->name = psStringCopy(name);
+        // This second copy of the name makes it easily available when writing the table
+        column->name = psStringCopy(name);
         switch (typecode) {
-           // TBYTE and TSHORT fall though to read into psS32
-          case TBYTE:
-          case TSHORT:
+            READ_TABLE_ROW_CASE(TBYTE, long, S32, S32);
+            READ_TABLE_ROW_CASE(TSHORT, short, S16, S16);
+            READ_TABLE_ROW_CASE(TUSHORT, unsigned short, U16, U16);
+            READ_TABLE_ROW_CASE(TULONG, unsigned long, U32, U32);
             READ_TABLE_ROW_CASE(TLONG, long, S32, S32);
             READ_TABLE_ROW_CASE(TLONGLONG, psS64, S64, S64);
@@ -216 +250 @@
               column->data.str = psAlloc(sizeof(psString) * table->numRows);
               column->type = PS_DATA_STRING;
+              column->elementSize = 0;
               for (int row = 0; row < table->numRows; row++) {
                   psString value = psStringAlloc(repeat);
@@ -227 +262 @@
                   }
                   column->data.str[row] = value;
+                  int len = strlen(value);
+                  if (len > column->elementSize) {
+                    column->elementSize = len;
+                  }
               }
               break;
@@ -245 +284 @@
 }
 
+// return the index for a given column
 int psFitsTableGetColumnIndex(psFitsTable *table, psString name)
 {
@@ -258 +298 @@
 }
 
+// return F32 value for given cell
 psF32 psFitsTableGetCellF32(psFitsTable *table, int row, int col) {
     return table->columns[col].data.F32[row];
 }
 
+// return value for a named column for a given row
 psF32 psFitsTableGetF32(bool *status, psFitsTable *table, psString name, int row) {
     bool mdok;
@@ -295 +337 @@
                     psFitsTableColumn *column = &table->columns[col];
                     switch (column->type) {
+                    case PS_DATA_S16:    
+                        column->data.S16[newRow] = column->data.S16[row];
+                        break;
+                    case PS_DATA_U16:    
+                        column->data.U16[newRow] = column->data.U16[row];
+                        break;
                     case PS_DATA_S32:    
                         column->data.S32[newRow] = column->data.S32[row];
+                        break;
+                    case PS_DATA_U32:    
+                        column->data.U32[newRow] = column->data.U32[row];
                         break;
                     case PS_DATA_S64:    
@@ -331 +382 @@
         }
     }
+
     // set the new number of rows
-    // XXX: should we do this? 
     table->numRows = newRow;
 
@@ -338 +389 @@
 }
 
+// Get the TFORM character, given a PS type
+static inline char getTForm(psDataType type)
+{
+    switch (type) {
+    case PS_TYPE_U8:
+    case PS_TYPE_S8:
+        return 'B';
+    case PS_TYPE_S16:
+        return 'I';
+    case PS_TYPE_S32:
+        return 'J';
+    case PS_TYPE_U64:
+    case PS_TYPE_S64:
+        return 'K';
+    case PS_TYPE_U16:
+        return 'U';
+    case PS_TYPE_U32:
+        return 'V';
+    case PS_TYPE_F32:
+        return 'E';
+    case PS_TYPE_F64:
+        return 'D';
+    case PS_TYPE_BOOL:
+        return 'L';
+    case PS_DATA_STRING:
+        return 'A';
+    default:
+        psError(PS_ERR_UNKNOWN, true, "Unknown type: %x\n", type);
+        return '?';
+    }
+}
+
+static bool fitsInsertTableNew(psFits* fits,          // FITS file
+                            const psMetadata* header, // FITS header to write
+                            psFitsTable *table,       // internal representation of the table
+                            const char *extname,      // Extension name to give table
+                            bool after,               // Write table after current extension?
+                            bool writeData            // Write data?
+                            )
+{
+    PS_ASSERT_FITS_NON_NULL(fits, false);
+    PS_ASSERT_FITS_WRITABLE(fits, false);
+    PS_ASSERT_PTR_NON_NULL(table, false);
+
+    int status = 0;
+
+    long numRows = table->numRows;
+    if (writeData && numRows < 1) {
+        // no table data, what can I do?
+        psError(PS_ERR_BAD_PARAMETER_SIZE, true,
+                _("Can't create a table without any rows."));
+        return false;
+    }
+
 #ifdef notdef
-
+    // Find the unique items in the array of metadata 'rows', and their sizes
+    psMetadata *colSpecs = psMetadataAlloc(); // Column specifications
+    size_t rowSize = 0;                 // Size (in bytes) of each row
+    for (long i = 0; i < numRows; i++) {
+        psMetadata* row = table->data[i];
+
+        if (!row) {
+            continue;
+        }
+        psMetadataIterator *rowIter = psMetadataIteratorAlloc(row, PS_LIST_HEAD, NULL); // Iterator
+        psMetadataItem *colItem;        // Column item, from iteration
+        while ((colItem = psMetadataGetAndIncrement(rowIter))) {
+            if (!(PS_DATA_IS_PRIMITIVE(colItem->type) || colItem->type == PS_DATA_STRING ||
+                    colItem->type == PS_DATA_VECTOR)) {
+                // unsupported type -- treating as an error
+                psError(PS_ERR_BAD_PARAMETER_TYPE, true,
+                        "Unsupported data type (%d) for Metadata Item '%s' in row %ld.",
+                        colItem->type, colItem->name, i);
+                psFree(rowIter);
+                psFree(colSpecs);
+                return false;
+            }
+
+            size_t size = columnSize(colItem); // Size for this column
+
+            // Check to see if we know about this one; or update the size if required
+            psMetadataItem *colSpecItem = psMetadataLookup(colSpecs, colItem->name);
+            if (!colSpecItem) {
+                // A new one!
+                colSpec *spec = psAlloc(sizeof(colSpec)); // Specification for this column
+                // BOOL type is not a valid vector type, so we translate it to U8
+                spec->type = colItem->type == PS_TYPE_BOOL ? PS_TYPE_U8 : colItem->type;
+                spec->size = size;
+                if (colItem->type == PS_DATA_VECTOR) {
+                    psVector *vector = colItem->data.V; // The vector
+                    spec->vectorType = vector->type.type;
+                }
+                psMetadataAddPtr(colSpecs, PS_LIST_TAIL, colItem->name, PS_DATA_UNKNOWN, "", spec);
+                psFree(spec);           // Drop reference
+                rowSize += PSELEMTYPE_SIZEOF(spec->type);
+            } else {
+                colSpec *spec = colSpecItem->data.V; // The specification
+                if (size > spec->size) {
+                    spec->size = size;
+                }
+                if (colItem->type != spec->type &&
+                    colItem->type != PS_TYPE_BOOL && spec->type != PS_TYPE_U8) {
+                    psWarning("Differing type found for column %s: %x vs %x --- using the first found.\n",
+                              colSpecItem->name, colItem->type, spec->type);
+                }
+                if (colItem->type == PS_DATA_VECTOR) {
+                    psVector *vector = colItem->data.V; // The vector
+                    if (vector->type.type != spec->vectorType) {
+                        psWarning("Differing vector type found for column %s: %x vs %x "
+                                 "--- using the first found.\n", colSpecItem->name, vector->type.type,
+                                 spec->vectorType);
+                    }
+                }
+            }
+        }
+        psFree(rowIter);
+    }
+
+    long numColumns = colSpecs->list->n;// Number of columns
+    if (numColumns == 0) {
+        // No table columns found
+        psError(PS_ERR_BAD_PARAMETER_SIZE, true,
+                "Did not find any column data to write to a table.");
+        psFree(colSpecs);
+        return false;
+    }
+#endif // notdef end of OLD column spec determination
+
+    // Create array of column names and types.
+    psArray *columnNames = psArrayAlloc(table->numCols); // Array of column names, for cfitsio
+    psArray *columnTypes = psArrayAlloc(table->numCols); // Array of column types, for cfitsio
+    for (long i = 0; i < table->numCols; i++) {
+        psFitsTableColumn *column = &table->columns[i];
+        columnNames->data[i] = psMemIncrRefCounter(column->name);
+        psString colType = NULL;        // The FITS column type
+        size_t size = column->elementSize;
+        if (column->type == PS_DATA_VECTOR) {
+            psStringAppend(&colType, "%zd%c", size, getTForm(column->vectorType));
+        } else {
+            psStringAppend(&colType, "%zd%c", size, getTForm(column->type));
+        }
+        columnTypes->data[i] = colType;
+    }
+
+    // Create the table HDU
+    int numHDUs = psFitsGetSize(fits);  // Number of HDUs in file
+    if (numHDUs == 0) {
+        // We're creating the first extension
+        fits_create_tbl(fits->fd,
+                        BINARY_TBL,
+                        writeData ? numRows : 0, // number of rows in table
+                        table->numCols, // number of columns in table
+                        (char**)columnNames->data, // names of the columns
+                        (char**)columnTypes->data, // format of the columns
+                        NULL, // physical unit of columns
+                        NULL, // skip extension name: we set the by hand below
+                        &status);
+    } else {
+        if (!after) {
+            if (psFitsGetExtNum(fits) == 0) {
+                // We're creating a replacement primary HDU.
+                // Set status to signal fits_insert_img to insert a new primary HDU
+                status = PREPEND_PRIMARY;
+            } else {
+                // Move back one to perform an insert after the previous HDU
+                psFitsMoveExtNum(fits, -1, true);
+            }
+        }
+        // Insert the table
+        fits_insert_btbl(fits->fd,
+                         writeData ? numRows : 0, // number of rows in table
+                         table->numCols, // number of columns in table
+                         (char**)columnNames->data, // names of the columns
+                         (char**)columnTypes->data, // format of the columns
+                         NULL, // physical unit of columns
+                         NULL, // skip extension name: we set this by hand below
+                         0, &status);
+    }
+    psFree(columnNames);
+    psFree(columnTypes);
+
+    if (status != 0) {
+        psFitsError(status, true, "Unable to create FITS table with %d columns and %d rows",
+                    table->numCols, table->numRows);
+        return false;
+    }
+
+    // Write header
+    if (header && !psFitsWriteHeader(fits, header)) {
+        psError(PS_ERR_IO, false, "Unable to write FITS header.\n");
+        return false;
+    }
+
+    // write the header, if any.
+    if (extname && strlen(extname) > 0) {
+        if (!psFitsSetExtName(fits, extname)) {
+            psError(PS_ERR_IO, false, "Unable to write FITS header extension name.\n");
+            return false;
+        }
+    }
+
+    if (writeData) {
+        // psMetadataIteratorSet(colSpecsIter, PS_LIST_HEAD);
+        for (long colNum = 1; colNum <= table->numCols; colNum++) {
+            // Note: colNum is unit-indexed, because it's for cfitsio
+            psFitsTableColumn *column = &table->columns[colNum - 1];
+            // colSpec *spec = colSpecItem->data.V; // The specification
+            if (PS_DATA_IS_PRIMITIVE(column->type)) {
+#ifdef notdef
+                size_t dataSize = PSELEMTYPE_SIZEOF(column->type); // Size (in bytes) of this type
+                psVector *columnData = psVectorAlloc(table->n, column->type); // The raw row data, to be written
+                psVectorInit(columnData, 0);
+                for (long i = 0; i < table->numRows; i++) {
+                    psMetadata *row = table->data[i]; // The row of interest
+                    psMetadataItem *dataItem = psMetadataLookup(row, colSpecItem->name); // Value of interest
+                    if (dataItem) {
+                        memcpy(&columnData->data.U8[i * dataSize], &dataItem->data, dataSize);
+                    } else {
+                        // this element is missing from this row; insert an appropriate-sized place holder
+                        // XXX this should insert a NAN for float / double and an appropriate blank for int types
+                        // XXX for the moment I am putting in 0.0
+                        memset(&columnData->data.U8[i * dataSize], 0, dataSize);
+                    }
+                }
+#endif
+
+                int fitsDataType;           // Data type for cfitsio
+                p_psFitsTypeToCfitsio(column->type, NULL, NULL, &fitsDataType);
+                fits_write_col(fits->fd,
+                               fitsDataType,
+                               colNum, // column number
+                               1, // first row
+                               1, // first element
+                               table->numRows, // number of rows
+                               column->data.U8, // the data
+                               &status);
+                // psFree(columnData);
+            } else {
+                switch (column->type) {
+                  case PS_DATA_STRING: {
+#ifdef notdef
+                      psArray *strings = psArrayAlloc(table->n); // Array of strings
+                      for (long i = 0; i < table->n; i++) {
+                          psMetadata *row = table->data[i]; // The row of interest
+                          strings->data[i] = psMemIncrRefCounter(psMetadataLookupStr(NULL, row,
+                                                                                     colSpecItem->name));
+                      }
+                      fits_write_col_str(fits->fd, colNum, 1, 1, table->n, (char**)strings->data, &status);
+#endif
+                      fits_write_col_str(fits->fd, colNum, 1, 1, table->numRows, (char**)column->data.str, &status);
+                      // psFree(strings);
+                      break;
+                  }
+                  case PS_DATA_VECTOR: {
+                      size_t dataSize = PSELEMTYPE_SIZEOF(column->vectorType); // Size of data, in bytes
+                      psVector *columnData = psVectorAlloc(column->elementSize * table->numRows * dataSize, PS_TYPE_U8);
+                      psVectorInit(columnData, 0);
+                      for (long i = 0; i < table->numRows; i++) {
+#ifdef notdef
+                          psMetadata *row = table->data[i]; // The row of interest
+                          psMetadataItem* dataItem = psMetadataLookup(row, colSpecItem->name);
+                          if (dataItem->type != PS_DATA_VECTOR) {
+                              // Just in case --- get a zero instead of some weird result
+                              continue;
+                          }
+                          psVector *vector = dataItem->data.V;
+#endif
+                          psVector *vector = column->data.vec[i];
+                          memcpy(&columnData->data.U8[i * dataSize * column->elementSize], vector->data.U8,
+                                 vector->n * dataSize);
+                      }
+
+                      int fitsDataType;           // Data type for cfitsio
+                      p_psFitsTypeToCfitsio(column->vectorType, NULL, NULL, &fitsDataType);
+                      fits_write_col(fits->fd, fitsDataType, colNum, 1, 1, table->numRows * column->elementSize,
+                                     columnData->data.U8, &status);
+                      psFree(columnData);
+                      break;
+                  }
+                  default:
+                    psAbort("Should never get here.\n");
+                }
+            }
+
+            // Check error status from writing column
+            if (status != 0) {
+                psFitsError(status, true, "Unable to write column %ld of FITS table", colNum);
+                return false;
+            }
+        }
+    }
+
+    // This forces a re-scan of the header to ensure everything's kosher.  We found this occassionally
+    // necessary for compressed images, which are tables, so perhaps it helps here too.  I guess it can't
+    // hurt.
+    ffrdef(fits->fd, &status);
+    if (psFitsError(status, true, "Could not re-scan HDU.")) {
+        return false;
+    }
+
+    return true;
+}
+
+// END OF new insert table
+
+bool psFitsWriteTableNew(psFits* fits,
+                      const psMetadata* header,
+                      psFitsTable* table,
+                      const char *extname)
+{
+    PS_ASSERT_FITS_NON_NULL(fits, false);
+    PS_ASSERT_FITS_WRITABLE(fits, false);
+    if (!psFitsMoveLast(fits)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to move to last extension to write table");
+        return false;
+    }
+    bool writeData = table->numRows > 0;
+    return fitsInsertTableNew(fits, header, table, extname, true, writeData);
+}
+
+#ifdef notdef
 
 psArray* psFitsReadTableColumn(const psFits* fits,

trunk/psLib/src/fits/psFitsTableNew.h

@@ -14 +14 @@
 
 typedef struct {
-    /* psString    name; */
+    psString    name;
     psDataType  type;
+    psDataType  vectorType;             // type inside if Vector
+    int         elementSize;            // 1 for primitives, strlen for strings
+                                        // vector length for vectors
     // only types supported by fits tables will be implemented
     union {
@@ -49 +52 @@
 
 psFitsTable *psFitsReadTableNew(const psFits *fits);
-bool psFitsWriteTableNew(psFits *fits, psFitsTable *table, const char *extname);
+bool psFitsWriteTableNew(psFits *fits, const psMetadata *header, psFitsTable *table, const char *extname);
 bool psFitsTableCensor(psFitsTable *table, bool *rowMask);