Changeset 28208 for trunk/psLib/src/fits/psFitsTable.c

Timestamp:

Jun 3, 2010, 5:18:11 PM (16 years ago)

Author:

Paul Price

Message:

Adding functions psFitsReadTableAllColumns and psFitsWriteTableAllColumns to more efficiently read/write a FITS table. Instead of a psMetadata for each row, the table is read as a psMetadata with each column implemented as a vector, array or image.

File:

• trunk/psLib/src/fits/psFitsTable.c (modified) (1 diff)

trunk/psLib/src/fits/psFitsTable.c

@@ -746 +746 @@
     return true;
 }
+
+
+psMetadata *psFitsReadTableAllColumns(const psFits *fits)
+{
+    PS_ASSERT_FITS_NON_NULL(fits, NULL);
+
+    if (!readTableCheck(fits)) {
+        return NULL;
+    }
+
+    int status = 0;                     // CFITSIO return status
+
+    long numRows = 0;                   // Number of rows in table
+    int numCols = 0;                    // Number of columns in table
+    fits_get_num_rows(fits->fd, &numRows, &status);
+    fits_get_num_cols(fits->fd, &numCols, &status);
+    if (psFitsError(status, true, "Failed to determine the size of the current HDU table.")) {
+        return NULL;
+    }
+
+    int hdutype;                        // Type of HDU: need to distinguish ASCII and binary tables
+    fits_get_hdu_type(fits->fd, &hdutype, &status);
+    if (psFitsError(status, true, "Could not determine the HDU type.")) {
+        return false;
+    }
+
+    psMetadata *table = psMetadataAlloc();     // Table to return
+    for (int col = 1; col <= numCols; col++) { // Fortran indexing
+        char name[FLEN_VALUE];           // Column name
+        if (hdutype == BINARY_TBL) {
+            fits_get_bcolparms(fits->fd, col, name,
+                               NULL, NULL, NULL, NULL, NULL, NULL, NULL, &status);
+        } else {
+            fits_get_acolparms(fits->fd, col, name,
+                               NULL, NULL, NULL, NULL, NULL, NULL, NULL, &status);
+        }
+
+        int cfitsioType = 0;             // Column type from CFITSIO
+        long repeat;                     // Number of repeats
+        fits_get_eqcoltype(fits->fd, col, &cfitsioType, &repeat, NULL, &status);
+        if (psFitsError(status, true, "Could not determine the column data for %s", name)) {
+            psFree(table);
+            return false;
+        }
+
+        psDataType pslibType = p_psFitsTypeFromCfitsio(cfitsioType); // Column type in psLib
+        if (pslibType == PS_DATA_STRING) {
+            // Strings
+            int width;                  // Width of strings
+            if (fits_get_col_display_width(fits->fd, col, &width, &status) != 0) {
+                psFitsError(status, true, "Could not determine the width of column %s", name);
+                psFree(table);
+                return NULL;
+            }
+            psArray *array = psArrayAlloc(numRows); // Array of strings from table
+            for (int i = 0; i < numRows; i++) {
+                array->data[i] = psStringAlloc(width);
+            }
+            fits_read_col_str(fits->fd, col, 1, 1, numRows, "", (char**)array->data, NULL, &status);
+            if (psFitsError(status, true, "Failed to read column %s", name)) {
+                psFree(array);
+                psFree(table);
+                return NULL;
+            }
+            if (!psMetadataAddArray(table, PS_LIST_TAIL, name, 0, NULL, array)) {
+                psError(PS_ERR_BAD_FITS, false, "Unable to add column %s to table", name);
+                psFree(array);
+                psFree(table);
+                return NULL;
+            }
+            psFree(array);
+        } else if (repeat == 1) {
+            // Single numbers
+            psVector* vector = psVectorAlloc(numRows, pslibType); // Vector from table
+            fits_read_col(fits->fd, cfitsioType, col, 1, 1, numRows, NULL,
+                          vector->data.U8, NULL, &status);
+            if (psFitsError(status, true, "Failed to read column %s", name)) {
+                psFree(vector);
+                psFree(table);
+                return NULL;
+            }
+            if (!psMetadataAddVector(table, PS_LIST_TAIL, name, 0, NULL, vector)) {
+                psError(PS_ERR_BAD_FITS, false, "Unable to add column %s to table", name);
+                psFree(vector);
+                psFree(table);
+                return NULL;
+            }
+            psFree(vector);
+        } else  {
+            // Vectors of numbers
+            psAssert(pslibType != PS_DATA_STRING, "Vectors of strings not handled");
+            psImage *image = psImageAlloc(repeat, numRows, pslibType); // Image to store vector of vectors
+            for (int row = 1, i = 0; row <= numRows; row++, i++) { // Fortran indexing for row
+                int anynul = 0;         // Any nulls in what was read?
+                fits_read_col(fits->fd, cfitsioType, col, row, 1, repeat, NULL,
+                              image->data.U8[i], &anynul, &status);
+                if (psFitsError(status, true, "Failed to read column %s row %d", name, row)) {
+                    psFree(image);
+                    psFree(table);
+                    return NULL;
+                }
+            }
+            if (!psMetadataAddImage(table, PS_LIST_TAIL, name, 0, NULL, image)) {
+                psError(PS_ERR_BAD_FITS, false, "Unable to add column %s to table", name);
+                psFree(image);
+                psFree(table);
+                return NULL;
+            }
+            psFree(image);
+        }
+    }
+
+    return table;
+}
+
+bool psFitsWriteTableAllColumns(psFits *fits, psMetadata *header, const psMetadata *table, const char *extname)
+{
+    PS_ASSERT_FITS_NON_NULL(fits, false);
+    PS_ASSERT_FITS_WRITABLE(fits, false);
+    PS_ASSERT_METADATA_NON_NULL(table, false);
+
+    psArray *columns = psListToArray(table->list); // Columns in psMetadataItems
+    int numCols = columns->n;                      // Number of columns
+    long numRows = 0;                              // Number of rows
+    psArray *names = psArrayAlloc(numCols);        // Column names
+    psArray *types = psArrayAlloc(numCols);        // Column types
+
+    for (int i = 0; i < numCols; i++) {
+        psMetadataItem *colItem = columns->data[i]; // Column
+        names->data[i] = psMemIncrRefCounter(colItem->name);
+
+        int size = 0;                   // Size of column
+        char tform = 0;                 // Type in character for FITS
+        long rows = 0;                   // Number of rows
+        switch (colItem->type) {
+          case PS_DATA_VECTOR:
+            size = 1;
+            psVector *vector = colItem->data.V; // Vector of interest
+            tform = getTForm(vector->type.type);
+            rows = vector->n;
+            break;
+          case PS_DATA_ARRAY:
+            tform = getTForm(PS_DATA_STRING);
+            psArray *array = colItem->data.V; // Array of interest
+            rows = array->n;
+            for (int i = 0; i < rows; i++) {
+                const char *string = array->data[i];
+                size = PS_MAX(size, strlen(string));
+            }
+            break;
+          case PS_DATA_IMAGE: ;
+            psImage *image = colItem->data.V; // Image of interest
+            tform = getTForm(image->type.type);
+            size = image->numCols;
+            rows = image->numRows;
+            break;
+          default:
+            psAbort("Unrecognised column type: %x", colItem->type);
+        }
+        psString type = NULL;           // Column type
+        psStringAppend(&type, "%zd%c", size, tform);
+        types->data[i] = type;
+        if (i == 0) {
+            numRows = rows;
+        } else if (numRows != rows) {
+            psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Column %d has differing size: %ld vs %ld",
+                    i, rows, numRows);
+            psFree(names);
+            psFree(types);
+            psFree(columns);
+            return false;
+        }
+    }
+
+    // Create the table HDU
+    int numHDUs = psFitsGetSize(fits);  // Number of HDUs in file
+    int status = 0;                     // Status from cfitsio
+    if (numHDUs == 0) {
+        // We're creating the first extension
+        fits_create_tbl(fits->fd, BINARY_TBL, numRows, numCols, (char**)names->data, (char**)types->data,
+                        NULL, NULL, &status);
+    } else {
+        fits_insert_btbl(fits->fd, numRows, numCols, (char**)names->data, (char**)types->data,
+                         NULL, NULL, 0, &status);
+    }
+    psFree(names);
+    psFree(types);
+    if (psFitsError(status, true, "Unable to create FITS table with %d columns and %ld rows",
+                    numCols, numRows)) {
+        psFree(columns);
+        return false;
+    }
+
+    // Write header
+    if (header && !psFitsWriteHeader(fits, header)) {
+        psError(psErrorCodeLast(), false, "Unable to write FITS header.\n");
+        psFree(columns);
+        return false;
+    }
+    if (extname && strlen(extname) > 0 && !psFitsSetExtName(fits, extname)) {
+        psError(psErrorCodeLast(), false, "Unable to write FITS header extension name.\n");
+        psFree(columns);
+        return false;
+    }
+
+    // cfitsio requires that we write the data by columns --- urgh!
+    for (long col = 1, i = 0; col <= numCols; col++, i++) {      // Fortran indexing
+        psMetadataItem *colItem = columns->data[i]; // Column
+        switch (colItem->type) {
+          case PS_DATA_VECTOR: {
+              psVector *vector = colItem->data.V; // Vector
+              int cfitsioType;                    // Data type for cfitsio
+              p_psFitsTypeToCfitsio(vector->type.type, NULL, NULL, &cfitsioType);
+              fits_write_col(fits->fd, cfitsioType, col, 1, 1, numRows, vector->data.U8, &status);
+              break;
+          }
+          case PS_DATA_ARRAY: {
+              psArray *array = colItem->data.V; // Array of strings
+              fits_write_col_str(fits->fd, col, 1, 1, numRows, (char**)array->data, &status);
+              break;
+          }
+          case PS_DATA_IMAGE: {
+              psImage *image = colItem->data.V;                                        // Image of interest
+              psDataType type = image->type.type;                                      // Type of data
+              psVector *vector = psVectorAlloc(image->numCols * image->numRows, type); // Vector from image
+              if (!p_psImageCopyToRawBuffer(vector->data.U8, image, type)) {
+                  psError(psErrorCodeLast(), false, "Unable to copy image to buffer");
+                  psFree(columns);
+                  return false;
+              }
+              int cfitsioType;            // Data type for cfitsio
+              p_psFitsTypeToCfitsio(type, NULL, NULL, &cfitsioType);
+              fits_write_col(fits->fd, cfitsioType, col, 1, 1, image->numRows * image->numCols,
+                             vector->data.U8, &status);
+              break;
+          }
+          default:
+            psAbort("Unrecognised column type: %x", colItem->type);
+        }
+        // Check error status from writing column
+        if (psFitsError(status, true, "Unable to write column %ld of FITS table", col)) {
+            psFree(columns);
+            return false;
+        }
+    }
+    psFree(columns);
+
+    // 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;
+}