Changeset 38329 for branches/eam_branches/ohana.20150429/src/libfits/matrix/F_compress_M.c

Timestamp:

May 28, 2015, 8:27:00 PM (11 years ago)

Author:

eugene

Message:

adding code for table compress/uncompress; image compress works (but is not quite compatible with cfitsio for some reason and byteswap is broken)

File:

• branches/eam_branches/ohana.20150429/src/libfits/matrix/F_compress_M.c (modified) (7 diffs)

branches/eam_branches/ohana.20150429/src/libfits/matrix/F_compress_M.c

@@ -3 +3 @@
 # include <zlib.h>
 
-# define ESCAPE { \
-  if (ztile != NULL) free (ztile); \
-  if (optname != NULL) { \
-    for (j = 0; j < Noptions; j++) { \
-      free (optname[j]); \
-      free (optvalue[j]); \
-    } \
-    free (optname); \
-    free (optvalue); \
-  } \
-  if (out != NULL) free (out); \
-  if (otile != NULL) free (otile); \
-  if (ntile != NULL) free (ntile); \
-  return (FALSE); }
-
-# define MOD_KEYWORD(ZNAME,NAME,TYPE,OUT,VALUE) { \
-    gfits_modify (theader, ZNAME, TYPE, 1, OUT);          \
-    gfits_delete (theader, NAME, 1); }
-
 // the user needs to specify the various compression options:
-// ztile[i] 
-// zcmptype
+// ztile[i] (may be NULL)
+// zcmptype 
 
 // steps to convert an image to a compressed table:
@@ -35 +16 @@
 // * update header keywords
 
-// XXX do not try to support compression of an image with an associated table (not valid)
-
-int gfits_compress_image (Header *header, Matrix *matrix, FTable *ftable, int *ztile, char *zcmptype) {
-
-  // off_t Nzdata, Nzrows, zcol;
-  // int i, j, status, zimage, Nout, max_tile_size;
-  // char cmptype[80];
-  // char zaxis[10], naxis[10], key[10], word[81], exttype[81], checksum[81], datasum[81];
-  // int Noptions, NOPTIONS, zblank, oblank;
-  // VarLengthColumn zdef;
-  // float zscale, zzero;
-
-  // int zdata_pixsize, odata_pixsize;
-  // int *otile = NULL;
-  // int *ntile = NULL;
-  // char **optname = NULL;
-  // char **optvalue = NULL;
-  // char *out = NULL;
-  // char *zdata = NULL;
+// we do not try to support compression of an image with an associated table (not valid)
+
+# define ESCAPE(A) { fprintf (stderr, "error in %s @ line %d\n", __func__, __LINE__); goto escape; }
+
+int gfits_compress_image (Header *header, Matrix *matrix, FTable *ftable, int *Ztile, char *zcmptype) {
+
+  int i;
+
+  // pointers to track the tile data
+  int *otile = NULL;
+  int *ntile = NULL;
+  char *raw = NULL;
+  char *zdata = NULL;
+
+  int *ztile = Ztile;
+
+  Header *theader = ftable->header;
 
   // determine the number of tiles (from ztile[] and image size)
@@ -67 +46 @@
 
   // creates an empty (Naxes = 2, Naxis[i] = 0) table header with XTENSION = BINTABLE, EXTNAME = COMPRESSED_IMAGE
-  gfits_create_table_header (&theader, "BINTABLE", "COMPRESSED_IMAGE");
+  if (!gfits_create_table_header (theader, "BINTABLE", "COMPRESSED_IMAGE")) ESCAPE(A);
 
   // by using "P" format, we are limited to 32bit pointers (2GB heap)
   // XXX how is the "B" format used?
-  gfits_define_bintable_column (&theader, "1PB(0)", "COMPRESSED_DATA", "compressed image data", "none", 1.0, 0.0);
-
-  // XXX I need to copy the bulk of the header from the matrix header -> ftable.header
-  gfits_create_table (&theader, ftable);
-  // copy original header to output header
-  // XXX gfits_copy_header (header, theader); -- I need to copy only the desired keywords.
+  if (!gfits_define_bintable_column (theader, "1PB(0)", "COMPRESSED_DATA", "compressed image data", "none", 1.0, 0.0)) ESCAPE (A);
+  if (!gfits_delete (theader, "TUNIT1", 1)) ESCAPE(A);
+
+  // allocates the default data array (ftable->buffer)
+  if (!gfits_create_table (theader, ftable)) ESCAPE (A);
+
+  // copy original header to output header (XXX this needs to be finished)
+  if (!gfits_copy_keywords_compress (header, theader)) ESCAPE (A);
   
-  // this allocates the data array 
+  // this allocates the pointer data array for the full set of tiles (NOT the heap)
   char *tmpbuffer = NULL;
   ALLOCATE_ZERO (tmpbuffer, char, Ntile);
-  gfits_set_bintable_column (&theader, &ftable, "COMPRESSED_DATA", tmpbuffer, Ntile);
+  if (!gfits_set_bintable_column (theader, ftable, "COMPRESSED_DATA", tmpbuffer, Ntile)) ESCAPE (A);
+  free (tmpbuffer);
 
   // add ZIMAGE from output header
-  gfits_modify_alt (theader, "ZIMAGE", "%t", 1, TRUE);
-
-  // translate image-specific keywords to Z versions (updates the Z version, removes the old version)
-  // XXX use the header or the structure as authoratative?
-  MOVE_KEYWORD ("ZBITPIX", "BITPIX", "%d");
-  MOVE_KEYWORD ("ZNAXIS",  "NAXIS",  "%d");
+  if (!gfits_modify_alt (theader, "ZIMAGE", "%t", 1, TRUE)) ESCAPE (A);
+
+  // define compression-specific keywords, update header as needed.
+  if (!gfits_modify (theader, "ZCMPTYPE", "%s", 1, zcmptype)) ESCAPE(A);
+
+  if (!gfits_modify_alt (theader, "ZSIMPLE", "%t", 1, TRUE)) ESCAPE (A);
+  // gfits_modify (header, "ZTENSION", "%s", 1, exttype);
+
+  // supply the Z* header values from the source image header
+  if (!gfits_modify (theader, "ZBITPIX", "%d", 1, header[0].bitpix)) ESCAPE (A);
+
+  if (!gfits_modify (theader, "ZNAXIS",  "%d", 1, header[0].Naxes)) ESCAPE (A);
+
+  char keyword[11];
+  for (i = 0; i < header[0].Naxes; i++) {
+    snprintf (keyword, 10, "ZNAXIS%d", i + 1);
+    if (!gfits_modify (theader, keyword, OFF_T_FMT, 1,  header[0].Naxis[i])) ESCAPE (A);
+  }
 
   for (i = 0; i < header->Naxes; i++) {
-    snprintf (zaxis, 10, "ZNAXIS%d", i + 1);
-    snprintf (naxis, 10, "NAXIS%d", i + 1);
-    MOD_KEYWORD (zaxis,  naxis,  OFF_T_FMT,  header->Naxis[i]);
-  }    
-
-  for (i = 0; i < header->Naxes; i++) {
-    snprintf (key, 10, "ZTILE%d", i + 1);
-    if (!gfits_modify (theader, key, "%d", 1, ztile[i])) ESCAPE;
-  }
-
-  // define compression-specific keywords, update header as needed.
-  if (!gfits_modify (header, "ZCMPTYPE", "%s", 1, zcmptype)) ESCAPE;
+    snprintf (keyword, 10, "ZTILE%d", i + 1);
+    if (!gfits_modify (theader, keyword, "%d", 1, ztile[i])) ESCAPE(A);
+  }
+
+  if (!gfits_modify (theader, "BSCALE", "%lf", 1, header[0].bscale)) ESCAPE (A);
+  if (!gfits_modify (theader, "BZERO",  "%lf", 1, header[0].bzero)) ESCAPE (A);
+
+  // other keywords to define:
+  // ZEXTEND
+  // ZBLOCKED
+  // ZHECKSUM
+  // ZDATASUM
+  // PCOUNT -> ZPCOUNT
+  // GCOUNT -> ZGCOUNT
+
+  float zscale = 1.0;
+  // gfits_modify (header, "ZSCALE", "%f", 1, zscale);
+
+  float zzero = 0.0;
+  // gfits_modify (header, "ZZERO", "%f", 1, zzero);
+
+  int zblank = 32767;
+  // gfits_modify (header, "ZBLANK", "%d", 1, zblank);
+
+  int oblank = 32767;
+  // gfits_modify (header, "BLANK", "%d", 1, oblank);
+
+# if (0)
+  // need to set ZNAMEnn and ZVALnn if they are defined and used
 
   // search for algorithm-specific keywords. these are used to control compression options
@@ -114 +125 @@
     if (!gfits_modify (header, key, "%s", 1, optvalue[i])) break;
   }
-
-# if (0)
-  // XXX get ZMASKCMP 
-  int zsimple, have_zsimple;
-  int have_ztension;
-
-  // NOTE target of %t must be int length
-  // XXX how do I specify this??
-  // have_zsimple  = gfits_scan_alt (header, "ZSIMPLE", "%t", 1, &zsimple);
-  // have_ztension = gfits_scan (header, "ZTENSION", "%s", 1, exttype);
-
-  // this is a very bogus case: we cannot have both keywords
-  if (have_zsimple && have_ztension) ESCAPE;
-
-  // if neither are present, we have an image that is not really following the standard:
-  // assume it is a PHU
-  if (!have_zsimple && !have_ztension) {
-    header->simple = TRUE;
-    gfits_extended_to_primary (header, header->simple, "Image data");
-
-    MOD_KEYWORD_ALT ("ZEXTEND",  "EXTEND",   "%t", &header->extend, header->extend);
-    MOD_KEYWORD_ALT ("ZBLOCKED", "BLOCKED",  "%t", &header->extend, header->extend);
-  }
-
-  // have_zsimple : image comes from a PHU
-  if (have_zsimple) {
-    header->simple = zsimple;
-    gfits_delete (header, "ZSIMPLE", 1);
-    gfits_extended_to_primary (header, header->simple, "Image data");
-
-    MOD_KEYWORD_ALT ("ZEXTEND",  "EXTEND",   "%t", &header->extend, header->extend);
-    MOD_KEYWORD_ALT ("ZBLOCKED", "BLOCKED",  "%t", &header->extend, header->extend);
-  } 
-
-  // have_ztension : image comes from an extension
-  if (have_ztension) {
-    gfits_delete (header, "ZTENSION", 1);
-    gfits_modify_extended (header, exttype, "Image extension");
-
-    // we may have an uncompressed PCOUNT / GCOUNT value, otherwise set to 0,1
-    if (gfits_scan (header, "ZPCOUNT", "%d", 1, &header->pcount)) {
-        gfits_delete (header, "ZPCOUNT", 1);
-        gfits_modify (header, "PCOUNT", "%d", 1, header->pcount);
-    } else {
-        header->pcount = 0;
-        gfits_modify (header, "PCOUNT", "%d", 1, header->pcount);
-    }
-    if (gfits_scan (header, "ZGCOUNT", "%d", 1, &header->gcount)) {
-        gfits_delete (header, "ZGCOUNT", 1);
-        gfits_modify (header, "GCOUNT", "%d", 1, header->gcount);
-    } else {
-        header->pcount = 1;
-        gfits_modify (header, "GCOUNT", "%d", 1, header->gcount);
-    }
-  } else {
-    header->pcount = 0;
-    header->gcount = 1;
-    gfits_modify (header, "PCOUNT", "%d", 1, header->pcount);
-    gfits_modify (header, "GCOUNT", "%d", 1, header->gcount);
-  }
-
-  MOD_KEYWORD ("ZHECKSUM", "CHECKSUM", "%s", checksum,        checksum);
-  MOD_KEYWORD ("ZDATASUM", "DATASUM",  "%s", datasum,         datasum);
-
-  zscale = 1;
-  gfits_scan (header, "ZSCALE", "%f", 1, &zscale);
-
-  zblank = 32767;
-  gfits_scan (header, "ZBLANK", "%d", 1, &zblank);
-
-  oblank = 32767;
-  gfits_scan (header, "BLANK", "%d", 1, &oblank);
-
-  zzero = 0;
-  gfits_scan (header, "ZZERO", "%f", 1, &zzero);
 # endif
 
-  // XXX how to define zdef?
   VarLengthColumn zdef;
-  if (!gfits_varlength_column_define (ftable, &zdef, zcol)) ESCAPE;
-
-  max_tile_size = gfits_imtile_maxsize (matrix, ztile);
-
-  // size of an element in the vartable heap section
-  zdata_pixsize = gfits_vartable_heap_pixsize (zdef.format);
-
-  // size of a pixel in the final image (not needed)
-  // raw_pixsize = gfits_uncompressed_data_pixsize (cmptype, header[0].bitpix, optname, optvalue, Noptions)
+  if (!gfits_varlength_column_define (ftable, &zdef, 1)) ESCAPE(A);
+
+  // allocate the working buffers to the max tile size
+  int max_tile_size = gfits_imtile_maxsize (matrix, ztile);
 
   // size of a pixel in the raw image tile:
-  raw_pixsize = abs(header[0].bitpix) / 8;
+  int raw_pixsize = abs(header[0].bitpix) / 8;
+
+  // size of a pixel in the raw image tile:
+  int tile_pixsize = gfits_uncompressed_data_pixsize (zcmptype, header[0].bitpix, NULL, NULL, 0);
 
   ALLOCATE (raw,   char, raw_pixsize*max_tile_size);
-  ALLOCATE (zdata, char, raw_pixsize*max_tile_size); // 
+  ALLOCATE (zdata, char, raw_pixsize*max_tile_size);
 
   // init the otile[] counters
-  gfits_imfile_start (matrix, otile);
+  if (!gfits_imtile_start (matrix, otile)) ESCAPE (A);
 
   // compress the data : copy into a tile, compress the tile, then add to the output table
@@ -217 +150 @@
 
     // size of the current tile
-    Nraw = gfits_tile_size (matrix, otile, ztile);
+    int Nraw = gfits_tile_size (matrix, otile, ztile);
 
     // copy the raw pixels from their native matrix locations to the temporary output buffer
-    if (!gfits_collect_data (matrix, raw, Nraw, raw_pixsize, otile, oblank, ztile, zblank, zscale, zzero)) return (FALSE);
+    if (!gfits_collect_data (matrix, raw, Nraw, raw_pixsize, otile, oblank, ztile, zblank, zscale, zzero)) ESCAPE(A);
 
     // XXX the tile must not be > 2GB
-
-    Nzdata = raw_pixsize*max_tile_size; // available space, replaced with actual output size on compression
-    if (!gfits_compress_data (zdata, &Nzdata, zcmptype, optname, optvalue, Noptions, raw, Nraw, raw_pixsize)) return (FALSE);
-
-    if (!gfits_varlength_column_add_data (ftable, zdata, Nzdata, i, &zdef)) return FALSE;
+   
+    // optname, optvalue = NULL, Noptions = 0
+
+    int Nzdata = raw_pixsize*max_tile_size; // available space, replaced with actual output size on compression
+    if (!gfits_compress_data (zdata, &Nzdata, zcmptype, NULL, NULL, 0, raw, Nraw, tile_pixsize)) ESCAPE(A);
+
+    if (!gfits_varlength_column_add_data (ftable, zdata, Nzdata, i, &zdef)) ESCAPE(A);
 
     // XXX need to activate
-    if (FALSE) gfits_byteswap_zdata (zdata, Nzdata, zdata_pixsize) return (FALSE);
+    // if (FALSE) gfits_byteswap_zdata (zdata, Nzdata, zdata_pixsize) ESCAPE(A);
 
     // update the otile[] counters, carrying to the next dimension if needed
-    gfits_imfile_next (matrix, ztile, ntile, otile);
-  }
+    gfits_imtile_next (matrix, ztile, ntile, otile);
+  }
+
+  if (!gfits_varlength_column_finish (ftable, &zdef)) ESCAPE(A);
+
+  free (raw);
+  free (zdata);
+
+  free (otile);
+  free (ntile);
+
+  if (!Ztile) free (ztile);
   return (TRUE);
+
+escape:
+  FREE (raw);
+  FREE (zdata);
+
+  FREE (otile);
+  FREE (ntile);
+
+  if (!Ztile) FREE (ztile);
+  return FALSE;
 }
 
-  // this takes place in three major steps:
-  // 1) rawdata -> tiledata (copy a tile chunk to the tiledata buffer)
-  // 2) compress tiledata -> zdata
-  // 3) add zdata to table data heap
-
-  // heapdata.pixsize : zdef.format
-  // zdata.pixsize : depends on compression
-  // odata.pixsize : depends on compression
-  // idata.pixsize : header.bitpix
-
+int gfits_copy_keywords_compress (Header *srchead, Header *tgthead) {
+
+  int i;
+
+  for (i = 0; i < srchead->datasize; i += FT_LINE_LENGTH) {
+  }
+  return TRUE;
+}
 
 // raw_pixsize is the bytes / pixel for the input matrix
@@ -298 +251 @@
       } } }
 
+  // we need to set up switches for all of the possible combinations:
+  // this macro is used at the inner switch to run the actual loop
+  // note we copy Optr to Iptr 
+# define SCALE_AND_DIST_FLOAT(TYPE, SIZE) {                                     \
+    TYPE *Optr = (TYPE *) &matrix->buffer[SIZE*(offset + start)];       \
+    for (j = 0; j < Ztile[0]; j++, Iptr++, Optr++) {                    \
+      if (!isfinite(*Optr)) {                                           \
+        *Iptr = zblank;                                                 \
+      } else {                                                          \
+        *Iptr = (*Optr - zzero) / zscale;                               \
+      } } }
+
   // this macro sets up the outer switch and calls above macro with all output bitpix options
 # define SETUP_INSIZE(TYPE, SIZE) {                     \
@@ -312 +277 @@
         break;                                          \
       case -32:                                         \
-        SCALE_AND_DIST (float, 4);                      \
+        SCALE_AND_DIST_FLOAT (float, 4);                        \
         break;                                          \
       case -64:                                         \
-        SCALE_AND_DIST (double, 8);                     \
+        SCALE_AND_DIST_FLOAT (double, 8);                       \
         break;                                          \
       default:                                          \