Changeset 15179 for trunk/psLib/src/fits/psFitsHeader.c

Timestamp:

Oct 3, 2007, 11:27:21 AM (19 years ago)

Author:

Paul Price

Message:

Merging branch with FITS compression development.

File:

• trunk/psLib/src/fits/psFitsHeader.c (modified) (15 diffs)

trunk/psLib/src/fits/psFitsHeader.c

@@ -7 +7 @@
  *  @author Robert DeSonia, MHPCC
  *
- *  @version $Revision: 1.34 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-08-10 02:23:42 $
+ *  @version $Revision: 1.35 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-10-03 21:27:21 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -19 +19 @@
 #include <unistd.h>
 
+#include "psAssert.h"
 #include "psFits.h"
 #include "string.h"
@@ -32 +33 @@
 
 #define MAX_STRING_LENGTH 256  // maximum length string for FITS routines
+#define NUM_EMPTY_KEYS 8                // Number of keywords before header is considered practically empty
 
 // list of FITS header keys to ignore; NULL-terminated
-static char* ignoreFitsKeys[] = {"", NULL};
+static const char* ignoreFitsKeys[] = { "", NULL};
 
 // List of FITS header keys that may be duplicated; NULL-terminated
-static char *duplicateFitsKeys[] = {"COMMENT", "HIERARCH", "HISTORY", NULL};
+static const char *duplicateFitsKeys[] = { "COMMENT", "HIERARCH", "HISTORY", NULL};
 
 // List of FITS header keys not to write (handled by cfitsio); NULL-terminated
-static char *noWriteFitsKeys[] = {"SIMPLE", "XTENSION", "BITPIX", "NAXIS", "EXTNAME", "BSCALE", "BZERO",
-                                  "TFIELDS", NULL};
+static const char *noWriteFitsKeys[] = { "SIMPLE", "XTENSION", "BITPIX", "NAXIS", "EXTNAME", "BSCALE",
+                                         "BZERO", "TFIELDS", "PCOUNT", "GCOUNT", "ZIMAGE", "ZBITPIX",
+                                         "ZCMPTYPE", NULL};
+
 // List of the start of FITS header keys not to write (handled by cfitsio); NULL-terminated
-static char *noWriteFitsKeyStarts[] = {"NAXIS", "TTYPE", "TFORM", NULL};
-
-// List of FITS header keys to be written with fits_write_comment
-
-psMetadata* psFitsReadHeader(psMetadata* out,
-                             const psFits* fits)
-{
-    if (fits == NULL) {
-        psError(PS_ERR_BAD_PARAMETER_NULL, true,
-                _("The input psFits object can not NULL."));
-        return NULL;
-    }
-
-    if (out == NULL) {
-        out = psMetadataAlloc();
-    }
+static const char *noWriteFitsKeyStarts[] = { "NAXIS", "TTYPE", "TFORM", "ZNAXIS", "ZTILE", "ZNAME", "ZVAL",
+                                              NULL};
+
+// List of FITS header keys that may be present if the header is considered "empty"; NULL-terminated
+static const char *emptyKeys[] = { "SIMPLE", "BITPIX", "NAXIS", "EXTEND", "COMMENT", "CHECKSUM", "DATASUM",
+                                   NULL };
+
+// Compare a keyword with a list of keywords; return true if it's in the list
+static bool keywordInList(const char *keyword, // Keyword to check
+                          const char **list // List of keywords
+    )
+{
+    for (const char **check = list; *check; ++check) {
+        if (strcmp(keyword, *check) == 0) {
+            return true;
+        }
+    }
+    return false;
+}
+
+
+bool psFitsCheckSingleCompressedImagePHU(const psFits *fits, psMetadata *header)
+{
+    PS_ASSERT_FITS_NON_NULL(fits, false);
+
+    if (!fits->conventions.compression) {
+        // User has turned off compression conventions; doesn't want any nasty surprises
+        return false;
+    }
+
+    if (psFitsGetExtNum(fits) != 0) {
+        // It's not the PHU, so it can't be the PHU for a single compressed image!
+        return false;
+    }
+
+    if (psFitsGetSize(fits) != 2) {
+        // No second extension, or multiple extensions
+        return false;
+    }
+
+    int numKeys;                        // Number of keywords in the header
+    int status = 0;                     // CFITSIO status
+    fits_get_hdrspace(fits->fd, &numKeys, 0, &status);
+    if (numKeys > NUM_EMPTY_KEYS) {
+        return false;
+    }
+
+    int bitpix, naxis;                  // Bits per pixel and number of axes
+    long naxes[MAX_COMPRESS_DIM];       // Dimensions
+    fits_get_img_param(fits->fd, MAX_COMPRESS_DIM, &bitpix, &naxis, naxes, &status);
+    if (naxis != 0) {
+        return false;
+    }
+
+    if (!header) {
+        for (int i = 1; i <= numKeys; i++) {
+            // Just want to read the keyword names, without parsing the values and stuffing into a metadata
+            char keyName[MAX_STRING_LENGTH];// Keyword name
+            char keyValue[MAX_STRING_LENGTH]; // Corresponding value
+            char keyComment[MAX_STRING_LENGTH]; // Corresponding comment
+            fits_read_keyn(fits->fd, i, keyName, keyValue, keyComment, &status);
+            if (!keywordInList(keyName, emptyKeys)) {
+                return false;
+            }
+        }
+    } else {
+        psMetadataIterator *iter = psMetadataIteratorAlloc(header, PS_LIST_HEAD, NULL); // Iterator
+        psMetadataItem *item;           // Item from iteration
+        while ((item = psMetadataGetAndIncrement(iter))) {
+            if (!keywordInList(item->name, emptyKeys)) {
+                psFree(iter);
+                return false;
+            }
+        }
+        psFree(iter);
+    }
+
+    if (!psFitsMoveExtNum(fits, 1, false)) {
+        psWarning("Unable to examine first extension as suspect compressed image.");
+        return false;
+    }
+
+    if (fits_is_compressed_image(fits->fd, &status)) {
+        return true;
+    }
+
+    // It's not a single compressed image PHU --- move back to the PHU for the user
+    if (!psFitsMoveExtNum(fits, 0, false)) {
+        psWarning("Unable to examine first extension as suspect compressed image.");
+        return false;
+    }
+
+    return false;
+}
+
+char *p_psFitsHeaderParseString(char *string)
+{
+    if (!string || strlen(string) == 0) {
+        return string;
+    }
+
+    char *fixed = string;       // Fixed version of the string
+    // remove the single-quotes at front/end
+    if (fixed[0] == '\'' && fixed[strlen(string)-1] == '\'') {
+        string[strlen(string)-1] = '\0'; // Remove the trailing quote
+        fixed += 1; // Advance past the leading quote
+    }
+    // Remove trailing spaces, which are not significant, according to the FITS standard
+    // http://archive.stsci.edu/fits/fits_standard/node31.html
+    char *lastSpace = NULL; // The last space in the string
+    while (strlen(fixed) > 1 && (lastSpace = strrchr(fixed, ' ')) && lastSpace[1] == '\0') {
+        // This is a trailing space, not a leading space.
+        lastSpace[0] = '\0'; // Truncate the string here
+    }
+
+    return fixed;
+}
+
+// Read the header
+static psMetadata *readHeader(const psFits *fits // FITS file from which to read header
+                              )
+{
+    assert(fits);
+
+    psMetadata *header = psMetadataAlloc(); // Header, to return
 
     // Get number of key names
@@ -74 +187 @@
 
         // Check to see if the keyword should be ignored
-        bool ignoreKey = false;         // Ignore this keyword?
-        for (int i = 0; ignoreFitsKeys[i] && !ignoreKey; i++) {
-            if (strcmp(keyName, ignoreFitsKeys[i]) == 0) {
-                ignoreKey = true;
-            }
-        }
-        if (ignoreKey) {
+        if (keywordInList(keyName, ignoreFitsKeys)) {
             // We're done here; skip to the next key
             continue;
         }
 
-#if 0
-        // This doesn't seem to be necessary
-        if (strncmp(keyName, "HIERARCH ", 9) == 0) {
-            char temp[MAX_STRING_LENGTH];
-            strcpy(temp, &keyName[9]);
-            strcpy(keyName, temp);
-        }
-#endif
-
         // Check to see if the keyword should be duplicated
         int dupFlag = 0;                // Duplicate flag
-        for (int i = 0; duplicateFitsKeys[i] && !dupFlag; i++) {
-            if (strcmp(keyName, duplicateFitsKeys[i]) == 0) {
-                dupFlag = PS_META_DUPLICATE_OK;
-            }
+        if (keywordInList(keyName, duplicateFitsKeys)) {
+            dupFlag = PS_META_DUPLICATE_OK;
         }
 
@@ -116 +212 @@
         case 'X': // bit
         case 'B': // byte
-            success = psMetadataAddS8(out, PS_LIST_TAIL, keyName, dupFlag, keyComment, atoi(keyValue));
+            success = psMetadataAddS8(header, PS_LIST_TAIL, keyName, dupFlag, keyComment, atoi(keyValue));
             break;
         case 'I': // short int.
@@ -122 +218 @@
             // Trap NAN, INF and -INF, which cfitsio doesn't handle.
             if (strncasecmp(keyValue, "NAN", 3) == 0) {
-                success = psMetadataAddF32(out, PS_LIST_TAIL, keyName, dupFlag, keyComment, NAN);
+                success = psMetadataAddF32(header, PS_LIST_TAIL, keyName, dupFlag, keyComment, NAN);
             } else if (strncasecmp(keyValue, "INF", 3) == 0) {
-                success = psMetadataAddF32(out, PS_LIST_TAIL, keyName, dupFlag, keyComment, INFINITY);
+                success = psMetadataAddF32(header, PS_LIST_TAIL, keyName, dupFlag, keyComment, INFINITY);
             } else if (strncasecmp(keyValue, "-INF", 4) == 0) {
-                success = psMetadataAddF32(out, PS_LIST_TAIL, keyName, dupFlag, keyComment, -INFINITY);
+                success = psMetadataAddF32(header, PS_LIST_TAIL, keyName, dupFlag, keyComment, -INFINITY);
             } else {
-                success = psMetadataAddS32(out, PS_LIST_TAIL, keyName, dupFlag, keyComment, atoi(keyValue));
+                success = psMetadataAddS32(header, PS_LIST_TAIL, keyName, dupFlag, keyComment,
+                                           atoi(keyValue));
             }
             break;
         case 'J': // int.
-            success = psMetadataAddS32(out, PS_LIST_TAIL, keyName, dupFlag, keyComment, atoi(keyValue));
+            success = psMetadataAddS32(header, PS_LIST_TAIL, keyName, dupFlag, keyComment, atoi(keyValue));
             break;
         case 'U': // unsigned int.
-            success = psMetadataAddU32(out, PS_LIST_TAIL, keyName, dupFlag, keyComment, atol(keyValue));
+            success = psMetadataAddU32(header, PS_LIST_TAIL, keyName, dupFlag, keyComment, atol(keyValue));
             break;
 
         case 'K': // long int. can't all fit in a psS32, put in psF64
         case 'F':
-            success = psMetadataAddF64(out, PS_LIST_TAIL, keyName, dupFlag, keyComment, atof(keyValue));
+            success = psMetadataAddF64(header, PS_LIST_TAIL, keyName, dupFlag, keyComment, atof(keyValue));
             break;
         case 'C': {
-                char *keyValueFixed = keyValue; // Fixed version of the string
-                // remove the single-quotes at front/end
-                if (keyValueFixed[0] == '\'' && keyValueFixed[strlen(keyValue)-1] == '\'') {
-                    keyValue[strlen(keyValue)-1] = '\0'; // Remove the trailing quote
-                    keyValueFixed += 1; // Advance past the leading quote
-                }
-                // Remove trailing spaces, which are not significant, according to the FITS standard
-                // http://archive.stsci.edu/fits/fits_standard/node31.html
-                char *lastSpace = NULL; // The last space in the string
-                while (strlen(keyValueFixed) > 1 && (lastSpace = strrchr(keyValueFixed, ' ')) &&
-                        lastSpace[1] == '\0') {
-                    // This is a trailing space, not a leading space.
-                    lastSpace[0] = '\0'; // Truncate the string here
-                }
-
-                // Need to trap NAN, INF and -INF written by psFitsWriteHeader.
-                // cfitsio won't write these, so we write them as strings, and then have to trap them on read.
-                if (strcasecmp(keyValueFixed, "NAN") == 0) {
-                    success = psMetadataAddF32(out, PS_LIST_TAIL, keyName, dupFlag, keyComment, NAN);
-                } else if (strcasecmp(keyValueFixed, "INF") == 0) {
-                    success = psMetadataAddF32(out, PS_LIST_TAIL, keyName, dupFlag, keyComment, INFINITY);
-                } else if (strcasecmp(keyValueFixed, "-INF") == 0) {
-                    success = psMetadataAddF32(out, PS_LIST_TAIL, keyName, dupFlag, keyComment, -INFINITY);
-                } else {
-                    success = psMetadataAddStr(out, PS_LIST_TAIL, keyName, dupFlag, keyComment,
-                                               keyValueFixed);
-                }
-                break;
+            char *keyValueFixed = p_psFitsHeaderParseString(keyValue); // Fixed version of the string
+
+            // Need to trap NAN, INF and -INF written by psFitsWriteHeader.
+            // cfitsio won't write these, so we write them as strings, and then have to trap them on read.
+            if (strcasecmp(keyValueFixed, "NAN") == 0) {
+                success = psMetadataAddF32(header, PS_LIST_TAIL, keyName, dupFlag, keyComment, NAN);
+            } else if (strcasecmp(keyValueFixed, "INF") == 0) {
+                success = psMetadataAddF32(header, PS_LIST_TAIL, keyName, dupFlag, keyComment, INFINITY);
+            } else if (strcasecmp(keyValueFixed, "-INF") == 0) {
+                success = psMetadataAddF32(header, PS_LIST_TAIL, keyName, dupFlag, keyComment, -INFINITY);
+            } else if (!fits->conventions.compression ||
+                       (strcmp(keyName, "EXTNAME") != 0 || strcmp(keyValueFixed, "COMPRESSED_IMAGE") != 0)) {
+                // Ignore EXTNAME=COMPRESSED_IMAGE if compression convention is to be respected
+                success = psMetadataAddStr(header, PS_LIST_TAIL, keyName, dupFlag, keyComment,
+                                           keyValueFixed);
             }
+            break;
+        }
         case 'L': {
-                bool temp = (keyValue[0] == 'T') ? 1 : 0;
-                success = psMetadataAddBool(out, PS_LIST_TAIL, keyName, dupFlag, keyComment, temp);
-                break;
-            }
+            bool temp = (keyValue[0] == 'T') ? 1 : 0;
+            success = psMetadataAddBool(header, PS_LIST_TAIL, keyName, dupFlag, keyComment, temp);
+            break;
+        }
         default:
             psError(PS_ERR_IO, true, _("Specified FITS metadata type, %c, is not supported."), keyType);
-            return out;
+            psFree(header);
+            return NULL;
         }
 
         if (!success) {
             psError(PS_ERR_UNKNOWN, false, _("Failed to add metadata item, %s."), keyName);
-            return out;
+            psFree(header);
+            return NULL;
         }
 
@@ -193 +281 @@
         (void)fits_get_errstatus(status, fitsErr);
         psError(PS_ERR_IO, true, _("Failed to add metadata item, %s."), fitsErr);
-        return false;
-    }
-
+        psFree(header);
+        return false;
+    }
+
+    return header;
+}
+
+
+psMetadata* psFitsReadHeader(psMetadata* out,
+                             const psFits* fits)
+{
+    PS_ASSERT_FITS_NON_NULL(fits, NULL);
+
+    psMetadata *header = readHeader(fits); // Header
+    if (!header) {
+        return NULL;
+    }
+
+    // Explore the potential case that this is an empty PHU, and the first extension contains the sole image,
+    // which is compressed.
+    if (psFitsCheckSingleCompressedImagePHU(fits, header)) {
+        // This is really what we want, not the empty PHU
+        psTrace("psLib.fits", 1,
+                "This PHU should really be a compressed image --- getting that header instead.");
+        psFree(header);
+        header = readHeader(fits);
+        if (!header) {
+            return NULL;
+        }
+    }
+
+    if (!out) {
+        return header;
+    }
+
+    // Need to move header onto the nominated metadata
+    psMetadataIterator *iter = psMetadataIteratorAlloc(header, PS_LIST_HEAD, NULL); // Iterator
+    psMetadataItem *item;           // Item from iteration
+    while ((item = psMetadataGetAndIncrement(iter))) {
+        // Need to look for MULTI, which won't be picked up using the iterator.
+        psMetadataItem *multiCheckItem = psMetadataLookup(header, item->name);
+        assert(multiCheckItem);
+        unsigned int flag = 0;      // Flag to indicate MULTI; otherwise default action
+        if (multiCheckItem->type == PS_DATA_METADATA_MULTI) {
+            flag = PS_META_DUPLICATE_OK;
+        }
+        if (!psMetadataAddItem(out, item, PS_LIST_TAIL, flag)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to add header item %s to extant metadata.",
+                    item->name);
+            psFree(iter);
+            psFree(header);
+            return NULL;
+        }
+    }
+    psFree(iter);
+    psFree(header);
     return out;
 }
@@ -201 +342 @@
 psMetadata* psFitsReadHeaderSet(psMetadata* out, const psFits* fits)
 {
-    if (fits == NULL) {
-        psError(PS_ERR_BAD_PARAMETER_NULL, true, _("The input psFits object can not NULL."));
-        psFree(out);
-        return NULL;
-    }
-
-    if (out == NULL) {
+    PS_ASSERT_FITS_NON_NULL(fits, NULL);
+
+    if (!out) {
         out = psMetadataAlloc();
     }
@@ -279 +416 @@
             // image, the NAXISn haven't been changed; or after converting to F32, the BITPIX hasn't been
             // changed) so we'll take care of that for them.
-            bool writeKey = true;           // Should we write this keyword?
-            for (int i = 0; noWriteFitsKeys[i] && writeKey; i++) {
-                if (strcmp(item->name, noWriteFitsKeys[i]) == 0) {
-                    writeKey = false;
-                }
-            }
-            if (!writeKey) {
+            if (keywordInList(item->name, noWriteFitsKeys)) {
                 // Don't write it; skip to the next key
                 continue;
@@ -294 +425 @@
                 // that go in are correct.  However, when we're writing a "blank" HDU (header only), we want
                 // to preserve NAXISn etc for reference, so we don't do this.
+                bool writeKey = true;   // Write this keyword?
                 for (int i = 0; noWriteFitsKeyStarts[i] && writeKey; i++) {
                     if (strncmp(item->name, noWriteFitsKeyStarts[i], strlen(noWriteFitsKeyStarts[i])) == 0) {
@@ -349 +481 @@
                         }
                     } else {
-                        fits_update_key(fits->fd, TFLOAT, item->name, &item->data.F32, item->comment, &status);
+                        fits_update_key(fits->fd, TFLOAT, item->name, &item->data.F32, item->comment,
+                                        &status);
                     }
                     break;
@@ -364 +497 @@
                         }
                     } else {
-                        fits_update_key(fits->fd, TDOUBLE, item->name, &item->data.F64, item->comment, &status);
+                        fits_update_key(fits->fd, TDOUBLE, item->name, &item->data.F64, item->comment,
+                                        &status);
                     }
                     break;
@@ -395 +529 @@
                       )
 {
-    if (!fits) {
-        psError(PS_ERR_BAD_PARAMETER_NULL, true, _("The input psFits object can not NULL."));
-        return false;
-    }
-
-    if (!output) {
-        psError(PS_ERR_BAD_PARAMETER_NULL, true, _("The input psMetadata was NULL.  Need a non-NULL psMetadata for operation to be performed."));
-        return false;
-    }
+    PS_ASSERT_FITS_NON_NULL(fits, false);
+    PS_ASSERT_METADATA_NON_NULL(output, false);
 
     return fitsWriteHeader(fits, output, true);
@@ -413 +540 @@
                      )
 {
-    if (!fits) {
-        psError(PS_ERR_BAD_PARAMETER_NULL, true, _("The input psFits object can not NULL."));
-        return false;
-    }
+    PS_ASSERT_FITS_NON_NULL(fits, false);
 
     // We allow output == NULL in order to write a minimal header.
@@ -458 +582 @@
 bool psFitsHeaderValidate(psMetadata *header)
 {
-    if (header == NULL) {
-        psError(PS_ERR_BAD_PARAMETER_NULL, true, _("The input psMetadata was NULL.  Need a non-NULL psMetadata for operation to be performed."));
-        return false;
-    }
+    PS_ASSERT_METADATA_NON_NULL(header, false);
 
     // Traverse the metadata list and inspect at each key