Changeset 15179 for trunk/psLib/src/fits/psFits.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/psFits.c (modified) (28 diffs)

trunk/psLib/src/fits/psFits.c

@@ -7 +7 @@
  *  @author Robert DeSonia, MHPCC
  *
- *  @version $Revision: 1.71 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-09-18 03:01:17 $
+ *  @version $Revision: 1.72 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-10-03 21:27:21 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -20 +20 @@
 
 #include "psFits.h"
+#include "psFitsHeader.h"
 #include "string.h"
 #include "psError.h"
@@ -65 +66 @@
             return false;
         }
-        fits->fd = NULL;
+        fits->fd = NULL;
     }
     return true;
@@ -74 +75 @@
     if (!fits) return;
     if (fits->fd) {
-        fitsClose(fits);
+        fitsClose(fits);
     }
     psFree (fits->extword);
@@ -81 +82 @@
 bool psFitsClose(psFits* fits)
 {
-    if (fits == NULL) {
-        psError(PS_ERR_BAD_PARAMETER_NULL, true,
-                _("The input psFits object can not NULL."));
-        return false;
-    }
+    PS_ASSERT_FITS_NON_NULL(fits, false);
 
     bool status = fitsClose(fits);
@@ -95 +92 @@
 psFits* psFitsOpen(const char* name, const char* mode)
 {
+    PS_ASSERT_STRING_NON_EMPTY(name, NULL);
+
     int status = 0;
     fitsfile *fptr = NULL;      /* Pointer to the FITS file */
-
-    if (name == NULL) {
-        psError(PS_ERR_BAD_PARAMETER_NULL, true,
-                _("Specified filename can not be NULL."));
-        return NULL;
-    }
 
     /* check the mode to determine how to open/create file */
@@ -131 +124 @@
         if (access(name, F_OK) == 0) {
             // file exists, delete old one first
-            remove
-                (name);
+            remove(name);
         }
 
@@ -140 +132 @@
         (void)fits_create_file
         #endif
-        (&fptr,
-         name,
-         &status);
+        (&fptr, name, &status);
         if (fptr == NULL || status != 0) {
-            char fitsErr[MAX_STRING_LENGTH]
-            ;
+            char fitsErr[MAX_STRING_LENGTH];
             fits_get_errstatus(status, fitsErr);
             psError(PS_ERR_IO, true,
@@ -158 +147 @@
         (void)fits_open_file
         #endif
-        (&fptr,
-         name,
-         iomode,
-         &status);
+        (&fptr, name, iomode, &status);
         if (fptr == NULL || status != 0) {
-            char fitsErr[MAX_STRING_LENGTH]
-            ;
+            char fitsErr[MAX_STRING_LENGTH];
             fits_get_errstatus(status, fitsErr);
             psError(PS_ERR_IO, true,
@@ -177 +162 @@
     fits->writable = (iomode == READWRITE);
     fits->extword = NULL;
+    fits->conventions.compression = true;
     psMemSetDeallocator(fits,(psFreeFunc)fitsFree);
 
@@ -182 +168 @@
 }
 
-static void psFitsOptionsFree(psFitsOptions *opt)
-{
-    psFree(opt->tilesize);
-}
-
-psFitsOptions* psFitsOptionsAlloc(
+psErrorCode p_psFitsError(const char* filename, unsigned int lineno, const char* func, int status,
+                          bool new, const char *errorMsg, ...)
+{
+    if (status == 0) {
+        return PS_ERR_NONE;
+    }
+
+    va_list ap;                         // Variable arguments
+    va_start(ap, errorMsg);
+    psString msg = NULL;                // Message to pass to psError
+    psStringAppendV(&msg, errorMsg, ap);
+    va_end(ap);
+
+    char cfitsioMsg[MAX_STRING_LENGTH];   // Error message from cfitsio
+    (void)fits_get_errstatus(status, cfitsioMsg);
+
+    psStringAppend(&msg, "[CFITSIO error: %s]", cfitsioMsg);
+
+    psErrorCode code = p_psError(filename, lineno, func, PS_ERR_IO, new, msg); // Error code
+    psFree(msg);
+    return code;
+}
+
+static void psFitsCompressionFree(psFitsCompression *comp)
+{
+    PS_ASSERT_PTR_NON_NULL(comp,);
+    psFree(comp->tilesize);
+}
+
+psFitsCompression* psFitsCompressionAlloc(
     psFitsCompressionType type,         ///< type of compression
     psVector *tilesize,                 ///< vector defining compression tile size
@@ -195 +205 @@
 )
 {
-    psFitsOptions *opt = psAlloc(sizeof(psFitsOptions));
-
-    opt->type = type;
-    opt->tilesize = psVectorCopy(NULL, tilesize, PS_DATA_S64);
-    opt->noisebits = noisebits;
-    opt->scale = scale;
-    opt->smooth = smooth;
-
-    psMemSetDeallocator(opt, (psFreeFunc) psFitsOptionsFree);
-
-    return opt;
+    psFitsCompression *comp = psAlloc(sizeof(psFitsCompression));
+    psMemSetDeallocator(comp, (psFreeFunc) psFitsCompressionFree);
+
+    comp->type = type;
+    comp->tilesize = psVectorCopy(NULL, tilesize, PS_DATA_S64);
+    comp->noisebits = noisebits;
+    comp->scale = scale;
+    comp->smooth = smooth;
+
+    return comp;
 }
 
@@ -224 +233 @@
 }
 
-// move to the first HDU where extword == extname.  this is equivalent to fits_movnam_hdu() for
-// a user-defined word in place of EXTNAME
-bool p_psFitsMoveExtName_UserKey(const psFits *fits,
-                                 const char *extname,
-                                 const char *extword)
-{
-    PS_ASSERT_PTR_NON_NULL(fits,    false);
-    PS_ASSERT_PTR_NON_NULL(extname, false);
-    PS_ASSERT_PTR_NON_NULL(extword, false);
-
-    int hdutype = 0;
-    char name[MAX_STRING_LENGTH];
-    char extstring[MAX_STRING_LENGTH];
-
-    sprintf (extstring, "'%s'", extname);
-
-    // NOTE: fits_* return 0 for success
-    for (int i = 1; true; i++) {
-        // are we able to read the next HDU?
-
-        int status = 0;
-        if (fits_movabs_hdu(fits->fd, i, &hdutype, &status)) {
-            char fitsErr[MAX_STRING_LENGTH];
-            fits_get_errstatus(status, fitsErr);
-            psError(PS_ERR_LOCATION_INVALID, true,
-                    _("Could not find HDU with %s = '%s'. CFITSIO Error: %s"),
-                    extword, extname, fitsErr);
-            return false;
-        }
-        // is there a keyword called 'extword'? (read as string regardless of type)
-        status = 0;
-        if (fits_read_keyword(fits->fd, (char *)extword, name, NULL, &status)) {
-            continue;
-        }
-        // if this was read as a string, we will have leading and trailing single-quotes
-        // try both for comparison
-        
-        // do we have the right hdu (names match)?
-        if (!strcmp (name, extname) || !strcmp (name, extstring)) {
-            return true;
-        }
-    }
-    psAbort("we should not reach here");
-}
-
-// XXX I will need to define a low-level function p_psFitsMoveExtName_UserKey () which
-// uses fits_movabs_hdu() to replicate the functionality of fits_movnam_hdu using an
-// alternate name for EXTNAME
+// Files compressed with cfitsio's "imcopy" program may have multiple EXTNAME keywords, with the first set to
+// COMPRESSED_IMAGE.  However, fits_movnam_hdu won't find the second (proper) value of EXTNAME, and so can
+// fail to find a perfectly legitimate extension, simply because imcopy does something silly.  However, we
+// really want to be able to read these files (MegaCam data are shipped as imcopy-compressed images).
+// Therefore, we implement our own version of moving to an extension specified by name.  The pure cfitsio
+// version is used if "conventions.compression" handling is turned off in the psFits structure.
 bool psFitsMoveExtName(const psFits* fits,
                        const char* extname)
 {
-    int status = 0;
-
-    if (fits == NULL) {
-        psError(PS_ERR_BAD_PARAMETER_NULL, true,
-                _("The input psFits object can not NULL."));
-        return false;
-    }
-
-    if (extname == NULL) {
-        psError(PS_ERR_BAD_PARAMETER_NULL, true,
-                _("Specified extension name can not be NULL."));
-        return false;
-    }
-
-    if (fits->extword != NULL) {
-        bool result = p_psFitsMoveExtName_UserKey(fits, extname, fits->extword);
-        return (result);
-    }
-
-    if (fits_movnam_hdu(fits->fd, ANY_HDU, (char*)extname, 0, &status) != 0) {
-        char fitsErr[MAX_STRING_LENGTH];
-        fits_get_errstatus(status, fitsErr);
-        psError(PS_ERR_LOCATION_INVALID, true,
-                _("Could not find HDU '%s'. CFITSIO Error: %s"),
-                extname, fitsErr);
-        return false;
-    }
-
-    return true;
+    PS_ASSERT_FITS_NON_NULL(fits, false);
+    PS_ASSERT_STRING_NON_EMPTY(extname, false);
+
+    int status = 0;
+
+    if (!fits->conventions.compression && !fits->extword) {
+        // User wants to use cfitsio.  Good luck to them!
+        if (fits_movnam_hdu(fits->fd, ANY_HDU, (char*)extname, 0, &status) != 0) {
+            psFitsError(status, true, _("Could not find HDU '%s'"), extname);
+            return false;
+        }
+        return true;
+    }
+
+    bool ignoreCI = (fits->conventions.compression &&
+                     (strcmp(extname, "COMPRESSED_IMAGE") != 0)); // Ignore COMPRESSED_IMAGE extension name?
+    char *extword = (fits->extword ? fits->extword : "EXTNAME"); // Word to use as extension name
+
+#if 0
+    // XXX Future optimisation: loop through from the current HDU to the end, then from the start to the
+    // current position.  This will save seeking through the file multiple times.
+    int currentExt = psFitsGetExtNum(fits); // Current extension number
+    int numExt = psFitsGetSize(fits);   // Total number of extensions
+#endif
+
+    for (int i = 1; true; i++) {
+        int hdutype = 0;
+        if (fits_movabs_hdu(fits->fd, i, &hdutype, &status)) {
+            // We've run off the end
+            psFitsError(status, true, _("Could not find HDU with %s = '%s'"), extword, extname);
+            return false;
+        }
+        // Is there a keyword called 'extword'? (read as string regardless of type)
+        char name[MAX_STRING_LENGTH];  // Name of extension
+        if (fits_read_keyword(fits->fd, extword, name, NULL, &status)) {
+            // It doesn't exist in the header.
+            // This isn't the extension you're looking for.  Move along.
+            status = 0;
+            continue;
+        }
+        char *fixed = p_psFitsHeaderParseString(name); // Parsed version (removing quotes and spaces)
+
+        if (ignoreCI && strcmp(fixed, "COMPRESSED_IMAGE") == 0) {
+            // Read it again, Sam
+            if (fits_read_keyword(fits->fd, extword, name, NULL, &status)) {
+                status = 0;
+                continue;
+            }
+            fixed = p_psFitsHeaderParseString(name);
+        }
+
+        if (strcmp(fixed, extname) == 0) {
+            // We've arrived
+            return true;
+        }
+    }
+    psAbort("Should never reach here.");
 }
 
@@ -310 +305 @@
                       bool relative)
 {
-    if (fits == NULL) {
-        psError(PS_ERR_BAD_PARAMETER_NULL, true,
-                _("The input psFits object can not NULL."));
-        return false;
-    }
+    PS_ASSERT_FITS_NON_NULL(fits, false);
 
     int status = 0;
@@ -346 +337 @@
 bool psFitsMoveLast(psFits* fits)
 {
-    if (fits == NULL) {
-        psError(PS_ERR_BAD_PARAMETER_NULL, true,
-                _("The input psFits object can not NULL."));
-        return false;
-    }
+    PS_ASSERT_FITS_NON_NULL(fits, false);
+
     int size = psFitsGetSize(fits);
     if (size == 0) { // empty file -- no action needed
@@ -361 +349 @@
 int psFitsGetExtNum(const psFits* fits)
 {
+    PS_ASSERT_FITS_NON_NULL(fits, false);
     int hdunum;
-
-    if (fits == NULL) {
-        psError(PS_ERR_BAD_PARAMETER_NULL, true,
-                _("The input psFits object can not NULL."));
-        return PS_FITS_TYPE_NONE;
-    }
-
-
     return fits_get_hdu_num(fits->fd,&hdunum) - 1;
 }
@@ -375 +356 @@
 psString psFitsGetExtName(const psFits* fits)
 {
-    if (fits == NULL) {
-        psError(PS_ERR_BAD_PARAMETER_NULL, true,
-                _("The input psFits object can not NULL."));
-        return NULL;
-    }
+    PS_ASSERT_FITS_NON_NULL(fits, NULL);
 
     int status = 0;
@@ -389 +366 @@
         psError(PS_ERR_BAD_PARAMETER_NULL, true,
                 _("Header keyword %s is not found"), extword);
-        return NULL;
+        return NULL;
     }
     return psStringCopy(name);
@@ -396 +373 @@
 bool psFitsSetExtName(psFits* fits, const char* name)
 {
-    if (fits == NULL) {
-        psError(PS_ERR_BAD_PARAMETER_NULL, true,
-                _("The input psFits object can not NULL."));
-        return false;
-    }
-
-    if (name == NULL) {
-        psError(PS_ERR_BAD_PARAMETER_NULL, true,
-                _("Specified extension name can not be NULL."));
-        return false;
-    }
+    PS_ASSERT_FITS_NON_NULL(fits, false);
+    PS_ASSERT_STRING_NON_EMPTY(name, false);
 
     int status = 0;
@@ -428 +396 @@
                         bool relative)
 {
-    if (fits == NULL) {
-        psError(PS_ERR_BAD_PARAMETER_NULL, true,
-                _("The input psFits object can not NULL."));
-        return false;
-    }
+    PS_ASSERT_FITS_NON_NULL(fits, false);
 
     if (! fits->writable) {
@@ -466 +430 @@
                          const char* extname)
 {
-    if (fits == NULL) {
-        psError(PS_ERR_BAD_PARAMETER_NULL, true,
-                _("The input psFits object can not NULL."));
-        return false;
-    }
+    PS_ASSERT_FITS_NON_NULL(fits, false);
+    PS_ASSERT_STRING_NON_EMPTY(extname, false);
 
     if (! fits->writable) {
@@ -504 +465 @@
 int psFitsGetSize(const psFits* fits)
 {
-    if (fits == NULL) {
-        psError(PS_ERR_BAD_PARAMETER_NULL, true,
-                _("The input psFits object can not NULL."));
-        return 0;
-    }
+    PS_ASSERT_FITS_NON_NULL(fits, 0);
 
     int num = 0;
@@ -527 +484 @@
 psFitsType psFitsGetExtType(const psFits* fits)
 {
-    if (fits == NULL) {
-        psError(PS_ERR_BAD_PARAMETER_NULL, true,
-                _("The input psFits object can not NULL."));
-        return PS_FITS_TYPE_NONE;
-    }
+    PS_ASSERT_FITS_NON_NULL(fits, PS_FITS_TYPE_NONE);
 
     int status = 0;
@@ -556 +509 @@
 bool psFitsTruncate(psFits* fits)
 {
-    if (fits == NULL) {
-        psError(PS_ERR_BAD_PARAMETER_NULL, true,
-                _("The input psFits object can not NULL."));
-        return PS_FITS_TYPE_NONE;
-    }
+    PS_ASSERT_FITS_NON_NULL(fits, NULL);
 
     if (! fits->writable) {
@@ -597 +546 @@
 )
 {
+    PS_ASSERT_FITS_NON_NULL(fits, false);
+
     // convert psFitsCompressionType to cfitsio compression types
     int comptype;
@@ -640 +591 @@
         psAbort("can't map (long) type to a psLib type");
     }
+    psFree(dim);
     // status check belongs to fits_set_tile_dim() call
     if (status) {
@@ -681 +633 @@
 }
 
-
-bool psFitsSetOptions(
+psFitsCompression *psFitsCompressionGet(psFits* fits)
+{
+    PS_ASSERT_FITS_NON_NULL(fits, NULL);
+
+    int status = 0;                     // cfitsio status
+
+    psFitsCompressionType type = psFitsCompressionGetType(fits);
+    if (type < 0) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to get compression type.");
+        return NULL;
+    }
+
+    psElemType tileType;                // Type corresponding to "long"
+    if (sizeof(long) == sizeof(psS64)) {
+        tileType = PS_TYPE_S64;
+    } else if (sizeof(long) == sizeof(psS32)) {
+        tileType = PS_TYPE_S32;
+    } else {
+        psAbort("can't map (long) type to a psLib type");
+    }
+
+    psVector *tiles = psVectorAlloc(3, tileType); // Tile sizes
+    if (fits_get_tile_dim(fits->fd, 3, (long*)tiles->data.U8, &status)) {
+        psFitsError(status, true, "Unable to get compression tile sizes.");
+        psFree(tiles);
+        return NULL;
+    }
+
+    int noisebits;                      // Noise bits for compression
+    if (fits_get_noise_bits(fits->fd, &noisebits, &status)) {
+        psFitsError(status, true, "Unable to get compression noise bits.");
+        psFree(tiles);
+        return NULL;
+    }
+
+    int hscale = 0, hsmooth = 0;        // Scaling and smoothing for HCOMPRESS
+
+#if FITS_HCOMP
+    if (fits_get_hcomp_scale(fits->fd, &hscale, &status)) {
+        psFitsError(status, true, "Unable to get HCOMPRESS scaling.");
+        psFree(tiles);
+        return NULL;
+    }
+    if (fits_get_hcomp_smooth(fits->fd, &hsmooth, &status)) {
+        psFitsError(status, true, "Unable to get HCOMPRESS smoothing.");
+        psFree(tiles);
+        return NULL;
+    }
+#endif // FITS_HCOMP
+
+    psFitsCompression *compress = psFitsCompressionAlloc(type, tiles, noisebits, hscale, hsmooth);
+    psFree(tiles);                      // Drop reference
+
+    return compress;
+}
+
+psFitsCompressionType psFitsCompressionGetType(psFits* fits)
+{
+    PS_ASSERT_FITS_NON_NULL(fits, -1);
+
+    int status = 0;                     // cfitsio status
+    int comptype = 0;                   // cfitsio compression type
+    if (fits_get_compression_type(fits->fd, &comptype, &status)) {
+        psFitsError(status, true, "Unable to get compression type.");
+        return -1;
+    }
+
+    psFitsCompressionType type;
+    switch (comptype) {
+      case 0:
+        type = PS_FITS_COMPRESS_NONE;
+        break;
+      case GZIP_1:
+        type = PS_FITS_COMPRESS_GZIP;
+        break;
+      case RICE_1:
+        type = PS_FITS_COMPRESS_RICE;
+        break;
+      case HCOMPRESS_1:
+        type = PS_FITS_COMPRESS_HCOMPRESS;
+        break;
+      case PLIO_1:
+        type = PS_FITS_COMPRESS_PLIO;
+        break;
+      default:
+        psError(PS_ERR_UNKNOWN, true, "cfitsio reports unknown compression type.");
+        return -1;
+    }
+
+    return type;
+}
+
+
+bool psFitsCompressionApply(
     psFits* fits,                       ///< psFits object to close
-    psFitsOptions *opt                  ///< options object
+    psFitsCompression *comp             ///< options object
 )
 {
-    return psFitsSetCompression(
-            fits,
-            opt->type,
-            opt->tilesize,
-            opt->noisebits,
-            opt->scale,
-            opt->smooth);
+    return psFitsSetCompression(fits, comp->type, comp->tilesize, comp->noisebits, comp->scale, comp->smooth);
 }
 
@@ -826 +864 @@
 
 
+psFitsCompressionType psFitsCompressionTypeFromString(const char *string)
+{
+    if (!string || strlen(string) == 0) {
+        psWarning("Unable to identify compression type --- none set.");
+        return PS_FITS_COMPRESS_NONE;
+    }
+
+    if (strcmp(string, "GZIP") == 0) return PS_FITS_COMPRESS_GZIP;
+    if (strcmp(string, "RICE") == 0) return PS_FITS_COMPRESS_RICE;
+    if (strcmp(string, "HCOMPRESS") == 0) return PS_FITS_COMPRESS_HCOMPRESS;
+    if (strcmp(string, "PLIO") == 0) return PS_FITS_COMPRESS_PLIO;
+
+    psWarning("Unable to identify compression type (%s) --- none set.", string);
+    return PS_FITS_COMPRESS_NONE;
+}