Changeset 38362 for branches/eam_branches/ohana.20150429/src/libfits/matrix/F_uncompress_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_uncompress_M.c (modified) (6 diffs)

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

@@ -3 +3 @@
 # include <zlib.h>
 # define VERBOSE_DUMP 0
+
+int gfits_distribute_gzp2 (Matrix *matrix, char *raw, int Nraw, int raw_bitpix, int *otile, int oblank, int *ztile, int zblank, float zscale, float zzero);
 
 # define ESCAPE { \
@@ -321 +323 @@
 
     // copy the uncompressed pixels into their correct locations            
-    if (!gfits_distribute_data (matrix, out, Nout, raw_bitpix, otile, oblank, ztile, zblank, zscale, zzero)) ESCAPE;
+    if (!strcasecmp(cmptype, "GZIP_2")) {
+      if (!gfits_distribute_gzp2 (matrix, out, Nout, raw_bitpix, otile, oblank, ztile, zblank, zscale, zzero)) ESCAPE;
+    } else {
+      if (!gfits_distribute_data (matrix, out, Nout, raw_bitpix, otile, oblank, ztile, zblank, zscale, zzero)) ESCAPE;
+    }
 
     if (VERBOSE_DUMP && (row == 0)) {
@@ -359 +365 @@
 int gfits_distribute_data (Matrix *matrix, char *raw, int Nraw, int raw_bitpix, int *otile, int oblank, int *ztile, int zblank, float zscale, float zzero) {
 
-  int i, j, start, offset, coord, Nline;
+  int i, j;
   int *counter = NULL;
   int *Ztile = NULL;
@@ -374 +380 @@
 
   // number of lines in the tile (in pixels)
-  Nline = 1;
+  int Nline = 1;
   for (i = 1; i < matrix->Naxes; i++) {
     Nline *= Ztile[i];
@@ -385 +391 @@
   // 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] + ...));
-  start = otile[matrix->Naxes-1]*ztile[matrix->Naxes-1];
+  int start = otile[matrix->Naxes-1]*ztile[matrix->Naxes-1];
   for (i = matrix->Naxes - 2; i >= 0; i--) {
-    coord = otile[i]*ztile[i];
+    int coord = otile[i]*ztile[i];
     start = start*matrix->Naxis[i] + coord;
   }
   
   // pixel offset in output array relative to tile start
-  offset = 0;
+  int offset = 0;
 
   // we need to set up switches for all of the possible combinations:
@@ -466 +472 @@
 }
 
+// bitpix is the input data size/type
+int gfits_distribute_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);
+  ALLOCATE (Ztile, int, matrix->Naxes);
+
+  // counter for current row in tile to copy
+  // 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];
+  }
+
+  // double check reported size (needs the Isize)
+  assert (Nraw == Ztile[0]*Nline);
+
+  // set the starting point of the tile:
+  // 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));
+        }
+        *srcptr = *rawptr; 
+        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);
+}
+