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

Timestamp:

Jun 3, 2015, 8:20:19 PM (11 years ago)

Author:

eugene

Message:

adding GZIP_2 and getting RICE_1 to work (mostly)

File:

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

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

@@ -3 +3 @@
 # include <zlib.h>
 # define VERBOSE_DUMP 0
+
+int gfits_collect_gzp2 (Matrix *matrix, char *raw, int Nraw, int raw_bitpix, int *otile, int oblank, int *ztile, int zblank, float zscale, float zzero);
 
 // the user needs to specify the various compression options:
@@ -73 +75 @@
   // define compression-specific keywords, update header as needed.
   if (!gfits_modify (theader, "ZCMPTYPE", "%s", 1, zcmptype)) ESCAPE;
+
+  // define compression-specific keywords, update header as needed.
+  if (!strcasecmp(zcmptype, "GZIP_1") || !strcasecmp(zcmptype, "GZIP_2")) {
+    if ((header[0].bitpix == -32) || (header[0].bitpix == -64)) {
+      if (!gfits_modify (theader, "ZQUANTIZ", "%s", 1, "NONE")) ESCAPE;
+    }
+  }
 
   if (!gfits_modify_alt (theader, "ZSIMPLE", "%t", 1, TRUE)) ESCAPE;
@@ -155 +164 @@
   int cmp_pixsize = gfits_compressed_data_pixsize (zcmptype, header[0].bitpix, NULL, NULL, 0);
 
+  // define compression-specific keywords, update header as needed.
+  if (!strcasecmp(zcmptype, "RICE_1")) {
+    if (!gfits_modify (theader, "ZNAME1", "%s", 1, "BLOCKSIZE")) ESCAPE;
+    if (!gfits_modify (theader, "ZVAL1", "%d", 1, 32)) ESCAPE;
+    if (!gfits_modify (theader, "ZNAME2", "%s", 1, "BYTEPIX")) ESCAPE;
+    if (!gfits_modify (theader, "ZVAL2", "%d", 1, raw_pixsize)) ESCAPE;
+  }
+
   // size of a pixel in the image tile -- this is probably not needed
   // int tile_pixsize = abs(header[0].bitpix) / 8;
@@ -176 +193 @@
     // copy the raw pixels from their native matrix locations to the temporary output buffer
     // for float -> int scaling by zscale, zzero may be applied
-    if (!gfits_collect_data (matrix, raw, Nraw, raw_bitpix, otile, oblank, ztile, zblank, zscale, zzero)) ESCAPE;
+    if (!strcasecmp (zcmptype, "GZIP_2")) {
+      if (!gfits_collect_gzp2 (matrix, raw, Nraw, raw_bitpix, otile, oblank, ztile, zblank, zscale, zzero)) ESCAPE;
+    } else {
+      if (!gfits_collect_data (matrix, raw, Nraw, raw_bitpix, otile, oblank, ztile, zblank, zscale, zzero)) ESCAPE;
+    }
 
     if (VERBOSE_DUMP && (i == 0)) {
@@ -214 +235 @@
     }
 
-    // all other compression modes require swapping after compression
-    // (compresssion & decompression operate on native ENDIAN)
-    if (strcasecmp(zcmptype, "GZIP_1")) {
+    // GZIP_1 is swapped before compression; GZIP_2 implies swapping.  All other
+    // compression modes require swapping after compression (compresssion & decompression
+    // operate on native ENDIAN)
+    if (strcasecmp(zcmptype, "GZIP_1") && strcasecmp(zcmptype, "GZIP_2")) {
       // Nzdata is number of bytes
       if (!gfits_byteswap_zdata (zdata, Nzdata, cmp_pixsize)) ESCAPE;
@@ -383 +405 @@
   return (TRUE);
 }
+
+// GZIP_2 organizes data in bytes by significance: (A1 A2 A3 A4), (B1 B2 B3 B4) -> (A1 B1), (A2, B2), (A3, B3), (A4, B4)
+// NOTE that this organization depends on the machine endian state. 
+// raw_pixsize is the bytes / pixel for the input matrix
+// place the raw image bytes for the current tile into the tile buffer
+// * otile is the counter for the current output tile
+// * ztile gives the size of the i-th dimension of the current tile
+// * raw_bitpix defines the size and type of the raw buffer
+int gfits_collect_gzp2 (Matrix *matrix, char *raw, int Nraw, int raw_bitpix, int *otile, int oblank, int *ztile, int zblank, float zscale, float zzero) {
+
+  int i, j, k;
+  int *counter = NULL;
+  int *Ztile = NULL;
+
+  ALLOCATE (counter, int, matrix->Naxes); // counter for current row in tile to copy
+  ALLOCATE (Ztile, int, matrix->Naxes); // true sizes of this tile (in pixels)
+  
+  for (i = 0; i < matrix->Naxes; i++) {
+    counter[i] = 0;
+    Ztile[i] = MIN ((matrix->Naxis[i] - otile[i]*ztile[i]), ztile[i]);
+  }
+
+  // number of lines in the tile (in pixels)
+  int Nline = 1;
+  int Npix = matrix->Naxis[0];
+  for (i = 1; i < matrix->Naxes; i++) {
+    Nline *= Ztile[i];
+    Npix  *= matrix->Naxis[i];
+  }
+
+  // confirm the tile size:
+  assert (Nraw == Ztile[0]*Nline);
+
+  // set the starting point of the tile in the matrix buffer
+  // start = otile[0]*ztile[0] + otile[1]*ztile[1]*Naxis[0] + otile[2]*ztile[2]*Naxis[0]*Naxis[1] + ...;
+  // start = otile[0]*ztile[0] + Naxis[0]*(otile[1]*ztile[1] + Naxis[1]*(otile[2]*ztile[2] + ...));
+  int start = otile[matrix->Naxes-1]*ztile[matrix->Naxes-1];
+  for (i = matrix->Naxes - 2; i >= 0; i--) {
+    int coord = otile[i]*ztile[i];
+    start = start*matrix->Naxis[i] + coord;
+  }
+  
+  int size = 0;
+  switch (raw_bitpix) {
+    case   8: size = 1; break;
+    case  16: size = 2; break;
+    case  32: size = 4; break;
+    case -32: size = 4; break;
+    case -64: size = 8; break;
+    default: abort();
+  }
+
+  // pixel offset in output array relative to tile start
+  int offset = 0;
+  
+  static int pass = 0;
+  for (k = 0; k < size; k++) {
+    for (i = 0; i < Nline; i++) {
+# ifdef BYTE_SWAP      
+      char *srcptr = &matrix->buffer[size*(offset + start) + (size - k - 1)];
+# else
+      char *srcptr = &matrix->buffer[size*(offset + start) + k];
+# endif
+      char *rawptr = &raw[i*size*Ztile[0] + k*Nraw];    
+      for (j = 0; j < Ztile[0]; j++, srcptr += size, rawptr ++) {               
+        if (FALSE && (i == 0) && (j < 4) && (pass == 0)) {
+          fprintf (stderr, "0x%02hhx -> 0x%02hhx\n", (int)(srcptr - matrix->buffer), (int)(rawptr - raw));
+        }
+        *rawptr = *srcptr; 
+        myAssert (srcptr - matrix->buffer < Npix*size, "memory overrun");
+        myAssert (rawptr - raw < Nraw*size, "memory overrun");
+      }
+
+      // update the counters, carrying to the next dimension if needed
+      for (j = 1; j < matrix->Naxes; j++) {
+        counter[j] ++;
+        if (counter[j] == Ztile[j]) {
+          counter[j] = 0;
+        } else {
+          break;
+        }
+      }
+      if (j == matrix->Naxes) assert (i == Nline - 1); // we should be done here...
+
+      // Naxes = 3
+      // offset = counter[1]*matrix->Naxis[0] + counter[2]*matrix->Naxis[0]*matrix->Naxis[1] + 
+      // offset = matrix->Naxis[0]*(counter[1] + matrix->Naxis[1]*(counter[2] + matrix->Naxis[2]*...))
+
+      // determine the offset of the next line relative to the start position
+      offset = counter[matrix->Naxes - 1];
+      for (j = matrix->Naxes - 2; j >= 0; j--) {
+        offset = offset*matrix->Naxis[j] + counter[j];
+      }      
+    }
+  }
+  pass ++;
+
+  free (counter);
+  free (Ztile);
+
+  return (TRUE);
+}