Changeset 41474 for trunk/Ohana/src/libfits/matrix/F_uncompress_M.c

Timestamp:

Jan 24, 2021, 3:37:52 PM (5 years ago)

Author:

eugene

Message:

clarify in comments the buffers and pixel sizes used in the different stages of compression/uncompression; catch unimplemented compression cases; fix error in pixel size for RICE1; cleanup verbose code; use correct sscanf format code in ASCII table reading code

File:

• trunk/Ohana/src/libfits/matrix/F_uncompress_M.c (modified) (12 diffs)

trunk/Ohana/src/libfits/matrix/F_uncompress_M.c

@@ -40 +40 @@
     gfits_modify (header, NAME, TYPE, 1, OUT); }
 
+# define DUMP_ROWBYTES(BUF, INFO)              \
+  if (VERBOSE_DUMP && (row == 0)) {            \
+    int k;                                     \
+    fprintf (stderr, "%s: ", INFO);            \
+    for (k = 0; k < 64; k++) {                 \
+      fprintf (stderr, "%02hhx", BUF[k]);    \
+      if (k % 4 == 3) fprintf (stderr, " ");   \
+    } fprintf (stderr, "\n"); }
+
 int gfits_uncompress_image (Header *header, Matrix *matrix, FTable *ftable) {
 
@@ -97 +106 @@
   } else {
     gfits_delete (header, "ZTILE1", 1);
+    // if ZTILE1 exists, all ZTILEn must exist:
     for (axis = 1; axis < header->Naxes; axis++) {
       snprintf (key, 10, "ZTILE%d", axis + 1); 
@@ -127 +137 @@
   }
 
-  // XXX get ZMASKCMP 
-  int zsimple, have_zsimple;
-  int have_ztension;
+  // check for ZMASKCMP (not yet supported 2021.01.23)
+  char zmaskcmp[81];
+  int have_zmaskcmp = gfits_scan (header, "ZMASKCMP", "%s", 1, zmaskcmp);
+  if (have_zmaskcmp) {
+    fprintf (stderr, "null pixels were compressed with %s, not implemented in Ohana\n", zmaskcmp);
+    ESCAPE;
+  }
 
   // NOTE target of %t must be int length
-  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
+  int zsimple;
+  int have_zsimple  = gfits_scan_alt (header, "ZSIMPLE", "%t", 1, &zsimple);
+  int have_ztension = gfits_scan (header, "ZTENSION", "%s", 1, exttype);
+
+  // this is a bogus case: we cannot have both keywords
   if (have_zsimple && have_ztension) ESCAPE;
 
@@ -179 +194 @@
     }
   } else {
+    // if this is a PHU (whether or not ZSIMPLE is present), PCOUNT & GCOUNT must be basic (0,1)
     header->pcount = 0;
     header->gcount = 1;
@@ -188 +204 @@
   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);
+  // get other basic keywords (default values supplied)
+  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);
 
   // find the COMPRESSED_DATA column (format should be 1PB, 1PI, 1PJ)
-  // is it required that this be the only column?
+  // XXX is it required that this be the only column? -- not if ZMASKCMP exists (not yet supported)
   int colnum;
   for (colnum = 1; TRUE; colnum++) {
@@ -242 +252 @@
 
   // heapdata.pixsize : zdef.format
-  // zdata.pixsize : depends on compression
-  // odata.pixsize : depends on compression
-  // idata.pixsize : header.bitpix
+  // zdata.pixsize : depends on compression method (cmp_pixsize below)
+  // odata.pixsize : depends on compression method (raw_pixsize below)
+  // idata.pixsize : depends on header.bitpix
 
   // size of an element in the vartable heap section
   // zdata_pixsize = gfits_vartable_heap_pixsize (zdef.format); XXX
 
-  // size of a pixel in the output from the decompression routine
+  // depending on the compression method, the data may need to be byteswapped either before
+  // or after uncompression:
+
+  // size of a pixel in the compressed data (cmp_pixsize), used for byteswap below
   int cmp_pixsize = gfits_compressed_data_pixsize (cmptype, header[0].bitpix, optname, optvalue, Noptions);
-  int raw_bitpix  = gfits_uncompressed_data_bitpix (cmptype, header[0].bitpix, optname, optvalue, Noptions);
-  int raw_pixsize = gfits_uncompressed_data_pixsize (cmptype, header[0].bitpix, optname, optvalue, Noptions);
+
+  // size of a pixel in the raw uncomprssed buffer (raw_pixsize) after uncompression (last option says RICE ignores bitpix)
+  int raw_pixsize = gfits_uncompressed_data_pixsize (cmptype, header[0].bitpix, optname, optvalue, Noptions, FALSE);
+
+  // we also need the bitpix value for the raw uncompressed buffer to distinguish int and floats
+  int raw_bitpix  = gfits_uncompressed_data_bitpix (cmptype, header[0].bitpix, optname, optvalue, Noptions, FALSE);
 
   // size of a pixel in the final image (not needed)
@@ -261 +278 @@
   ALLOCATE (out, char, Nout_alloc);
 
-  off_t row;
-
-  if (VERBOSE_DUMP) { 
-    int k;
-    fprintf (stderr, "unc tbl: "); 
-    for (k = 0; k < 64; k++) { fprintf (stderr, "%02hhx", ftable->buffer[k]); if (k % 4 == 3) fprintf (stderr, " "); } 
-    fprintf (stderr, "\n");
-  }
-
   // uncompress the data
   Nzrows = ftable->header->Naxis[1];
-  for (row = 0; row < Nzrows; row++) {
+  for (off_t row = 0; row < Nzrows; row++) {
+
+    DUMP_ROWBYTES (ftable->buffer, "tbl dat");
 
     // expected output size for this tile
@@ -280 +290 @@
     if (!zdata) ESCAPE;
 
-    if (VERBOSE_DUMP && (row == 0)) {
-      int k;
-      fprintf (stderr, "unc dat: "); 
-      for (k = 0; k < 64; k++) { fprintf (stderr, "%02hhx", zdata[k]); if (k % 4 == 3) fprintf (stderr, " "); } 
-      fprintf (stderr, "\n");
-    }
-
-    if (strcasecmp(cmptype, "NONE_2") && // strcasecmp(cmptype, "NONE") && 
-        strcasecmp(cmptype, "GZIP_1") && 
-        strcasecmp(cmptype, "GZIP_2") && 
-        strcasecmp(cmptype, "RICE_1") && 
-        strcasecmp(cmptype, "RICE_ONE")) {
-      // Nzdata is number of bytes
-      if (!gfits_byteswap_zdata (zdata, Nzdata, cmp_pixsize)) ESCAPE;
-    }
-
-    if (VERBOSE_DUMP && (row == 0)) {
-      int k;
-      fprintf (stderr, "unc SWP: "); 
-      for (k = 0; k < 64; k++) { fprintf (stderr, "%02hhx", zdata[k]); if (k % 4 == 3) fprintf (stderr, " "); } 
-      fprintf (stderr, "\n");
-      fprintf (stderr, "Nout  : %lu -> ", Nout);
-    }
+    DUMP_ROWBYTES (zdata, "cmp dat");
+
+    // byteswap the compressed data before uncompressing.
+    // compression algorithms which do not need to byteswap here return 1 for cmp_pixsize
+    if (!gfits_byteswap_zdata (zdata, Nzdata, cmp_pixsize)) ESCAPE;
+
+    DUMP_ROWBYTES (zdata, "cmp swp");
     
     // gfits_uncompress_data uncompresses from zdata to the temporary output buffer which must be allocated
-    // XXX the tile must not be > 2GB
-    // Nout going in is number of bytes 
+    // Note: the tile must not be > 2GB
+    // Nout going in is the expected number of bytes 
     if (!gfits_uncompress_data ((char *)zdata, Nzdata, cmptype, optname, optvalue, Noptions, out, &Nout, Nout_alloc, raw_pixsize)) ESCAPE;
 
-    if (VERBOSE_DUMP && (row == 0)) {
-      int k;
-      fprintf (stderr, "%lu\n", Nout);
-      fprintf (stderr, "unc swp: "); 
-      for (k = 0; k < 64; k++) { fprintf (stderr, "%02hhx", out[k]); if (k % 4 == 3) fprintf (stderr, " "); } 
-      fprintf (stderr, "\n");
-    }
+    DUMP_ROWBYTES (out, "unc swp");
     
     if (!strcasecmp(cmptype, "GZIP_1")) {
@@ -322 +310 @@
     }
     
-    if (VERBOSE_DUMP && (row == 0)) {
-      int k;
-      fprintf (stderr, "unc raw: "); 
-      for (k = 0; k < 64; k++) { fprintf (stderr, "%02hhx", out[k]); if (k % 4 == 3) fprintf (stderr, " "); } 
-      fprintf (stderr, "\n");
-    }
+    DUMP_ROWBYTES (out, "unc raw");
 
     // copy the uncompressed pixels into their correct locations            
@@ -336 +319 @@
     }
 
-    if (VERBOSE_DUMP && (row == 0)) {
-      int k;
-      fprintf (stderr, "unc mat: "); 
-      for (k = 0; k < 64; k++) { fprintf (stderr, "%02hhx", matrix->buffer[k]); if (k % 4 == 3) fprintf (stderr, " "); } 
-      fprintf (stderr, "\n");
-    }
+    DUMP_ROWBYTES (matrix->buffer, "out raw");
 
     // update the tile counters, carrying to the next dimension if needed
@@ -370 +348 @@
 }
 
-// bitpix is the input data size/type
+// raw_bitpix is the data size/type for the temporary buffer of uncompressed data
+// matrix->bitpix is the data size/type of the output image data
 int gfits_distribute_data (Matrix *matrix, char *raw, unsigned long int Nraw, int raw_bitpix, int *otile, int oblank, unsigned long int *ztile, int zblank, float zscale, float zzero) {
 
@@ -514 +493 @@
 }
 
-// bitpix is the input data size/type
+// raw_bitpix is the data size/type for the temporary buffer of uncompressed data
+// matrix->bitpix is the data size/type of the output image data
+// The gzip2 algorithm organizes the compressed data in groups by significant byte sequence (most significant byte of all pixels first, etc).
+// gfits_distribute_gzp2 needs to shuffle the bytes back into the correct location in the tile.
 int gfits_distribute_gzp2 (Matrix *matrix, char *raw, unsigned long int Nraw, int raw_bitpix, int *otile, int oblank, unsigned long int *ztile, int zblank, float zscale, float zzero) {
+  // unused parameters are supplied to make this function have the same API as gfits_distribute_data()
   OHANA_UNUSED_PARAM(oblank);
   OHANA_UNUSED_PARAM(zblank);