Changeset 15439

Timestamp:

Nov 2, 2007, 3:44:27 PM (19 years ago)

Author:

eugene

Message:

adding code from CFITSIO for reading various compressions

Location:

branches/eam_branch_20071015/Ohana/src/libfits

Files:

• Makefile (modified) (4 diffs)
• extern (added)
• extern/fits_hcompress.c (added)
• extern/fits_hdecompress.c (added)
• extern/pliocomp.c (added)
• extern/ricecomp.c (added)
• extern/ricecomp.h (added)
• header/F_copy_H.c (modified) (1 diff)
• include/gfitsio.h (modified) (2 diffs)
• matrix/F_compress_M.c (modified) (7 diffs)
• matrix/F_uncompress_data.c (added)
• table/F_table_varlength.c (modified) (6 diffs)

branches/eam_branch_20071015/Ohana/src/libfits/Makefile

@@ -10 +10 @@
 MATR    =       $(HOME)/matrix
 TABL    =       $(HOME)/table
+EXT     =       $(HOME)/extern
 INC     =       $(HOME)/include
 MAN     =       $(HOME)/doc
@@ -16 +17 @@
 # programs may add their own internal requirements here
 FULL_CFLAGS   = $(BASE_CFLAGS) -fPIC -Wall -Werror
-FULL_CPPFLAGS = $(BASE_CPPFLAGS)
+FULL_CPPFLAGS = $(BASE_CPPFLAGS) -I$(EXT)
 FULL_LDFLAGS  = $(BASE_LDFLAGS)
 
@@ -43 +44 @@
 $(MATR)/F_convert_format.$(ARCH).o              $(MATR)/F_read_segment.$(ARCH).o \
 $(MATR)/F_read_portion.$(ARCH).o                $(MATR)/F_load_M.$(ARCH).o       \
-$(MATR)/F_matrix.$(ARCH).o                      $(MATR)/F_compress_M.$(ARCH).o
+$(MATR)/F_matrix.$(ARCH).o                      $(MATR)/F_compress_M.$(ARCH).o   \
+$(MATR)/F_uncompress_data.$(ARCH).o
 
 TABLE_OBJ = \
@@ -54 +56 @@
 $(TABL)/F_table_varlength.$(ARCH).o
 
-OBJS = $(HEADER_OBJ) $(MATRIX_OBJ) $(TABLE_OBJ)
+EXTERN_OBJ = \
+$(EXT)/fits_hcompress.$(ARCH).o \
+$(EXT)/fits_hdecompress.$(ARCH).o \
+$(EXT)/pliocomp.$(ARCH).o \
+$(EXT)/ricecomp.$(ARCH).o
+
+OBJS = $(HEADER_OBJ) $(MATRIX_OBJ) $(TABLE_OBJ) $(EXTERN_OBJ)
+
+extern: $(EXTERN_OBJ)
 
 $(OBJS): $(INCS)

branches/eam_branch_20071015/Ohana/src/libfits/header/F_copy_H.c

@@ -8 +8 @@
 
   out[0].simple = in[0].simple;
+  out[0].unsign = in[0].unsign;
+  out[0].extend = in[0].extend;
   out[0].bitpix = in[0].bitpix;
-  out[0].extend = in[0].extend;
-  out[0].unsign = in[0].unsign;
-  out[0].size   = in[0].size;
-  out[0].bscale = in[0].bscale;
-  out[0].bzero  = in[0].bzero;
+
   out[0].Naxes  = in[0].Naxes;
   for (i = 0; i < FT_MAX_NAXES; i++) 
     out[0].Naxis[i] = in[0].Naxis[i];
+
+  out[0].size   = in[0].size;
+
+  out[0].pcount = in[0].pcount;
+  out[0].gcount = in[0].gcount;
+  out[0].bzero  = in[0].bzero;
+  out[0].bscale = in[0].bscale;
 
   ALLOCATE (out[0].buffer, char, out[0].size);

branches/eam_branch_20071015/Ohana/src/libfits/include/gfitsio.h

@@ -89 +89 @@
 
 typedef struct {
-  int  length;
-  int  offset;
-  int  maxlen;
-  int  nbytes;
-  int  Nstart;
-  int  theap;
-  char format;
-  char cmptype[81];
+  int  maxlen;                // max size of all table rows
+  int  nbytes;                // number of bytes per column element
+  int  Nstart;                // byte offset of this column
+  int  heap_start;            // byte offset to start of HEAP
+  char format;                // data format character (one of: XLABIJEDCM)
 } VarLengthColumn;
 
@@ -156 +153 @@
 int     gfits_write_matrix             PROTO((char *filename, Matrix *matrix)); 
 int     gfits_uncompress_image         PROTO((Header *header, Matrix *matrix, FTable *ftable, int primary));
-int     gfits_uncompress_data          PROTO((char *zdata, int Nzdata, char *cmptype, char **optname, char **optvalue, int Nopt, char *outdata, int *Nout));
+int     gfits_uncompress_data          PROTO((char *zdata, int Nzdata, char *cmptype, char **optname, char **optvalue, int Nopt, char *outdata, int *Nout, int out_pixsize));
 int     gfits_distribute_data          PROTO((Matrix *matrix, int bitpix, char *data, int Ndata, int *otile, int *ztile, float zscale, float zzero));
 int     gfits_byteswap_zdata           PROTO((char *zdata, int Nzdata, int bitpix));
+int     gfits_extension_is_compressed  PROTO((Header *header));
+int     gfits_tile_size                PROTO((Matrix *matrix, int *otile, int *ztile));
 
 /******************************* Table functions *************/

branches/eam_branch_20071015/Ohana/src/libfits/matrix/F_compress_M.c

@@ -8 +8 @@
   gfits_print (header, NAME, TYPE, 1, OUT); }
 
-# define MOD_KEYWORD_REQUIRED(NAME,ZNAME,TYPE,IN,OUT) { \
+# define MOD_KEYWORD_REQUIRED(ZNAME,NAME,TYPE,IN,OUT) { \
   if (!gfits_scan (header, ZNAME, TYPE, 1, IN)) ESCAPE; \
   gfits_delete (header, ZNAME, 1); \
@@ -30 +30 @@
 int gfits_uncompress_image (Header *header, Matrix *matrix, FTable *ftable, int primary) {
 
-  int i, j, status, zimage, zcol, Nzrows, Nout, Nzdata;
+  int i, j, status, zimage, zcol, Nzrows, Nout, Nzdata, max_tile_size, bytes_per_pixel;
   char cmptype[80];
   char zaxis[10], naxis[10], key[10], word[81], exttype[81], checksum[81], datasum[81];
@@ -145 +145 @@
   gfits_create_matrix (header, matrix);
 
-  // XXX is maxlen in bytes or in pixels?
-  ALLOCATE (out, char, zdef.maxlen);
-
   // counters for tile numbers
   ALLOCATE (otile, int, matrix->Naxes);
   ALLOCATE (ntile, int, matrix->Naxes);
+  max_tile_size = 1;
   for (i = 0; i < matrix->Naxes; i++) {
     otile[i] = 0;
     ntile[i] = (matrix->Naxis[i] % ztile[i]) ? (matrix->Naxis[i] / ztile[i] + 1) : (matrix->Naxis[i] / ztile[i]);
-  }
+    max_tile_size *= ztile[i];
+  }
+
+  // max_tile_size is in pixels
+  bytes_per_pixel = abs(header[0].bitpix) / 8;
+  ALLOCATE (out, char, bytes_per_pixel*max_tile_size);
 
   // uncompress the data
   Nzrows = ftable->header->Naxis[1];
   for (i = 0; i < Nzrows; i++) {
+
+    // expected output size for this tile
+    Nout = gfits_tile_size (matrix, otile, ztile);
+
     switch (zdef.format) {
-
-      // may need to multiple by byte/pix
-      Nout = zdef.maxlen;
-
       // gfits_varlength returns a pointer to a chunk of Nzdata of data elements starting at zdata.
       // uncompress_data uncompresses from zdata to the temporary output buffer which must be allocated
+
       // XXX need to have an API to send the user data
+
+      // XXX need to byte-swap the table column; this needs to be worked out more clearly
+      // in the APIs: do the table read / column extract functions swap or not?.  we have
+      // put the byte-swapping in gifts_varlength_column_pointer, but this is fairly weak.
 
       case 'B': {
         char *zdata;
         zdata = gfits_varlength_column_pointer (ftable, &zdef, i, &Nzdata);
-        gfits_byteswap_zdata (zdata, Nzdata, 8);
-        gfits_uncompress_data (zdata, Nzdata, zdef.cmptype, optname, optvalue, Noptions, out, &Nout);
-        gfits_distribute_data (matrix, 8, out, Nout, otile, ztile, zscale, zzero);      // copy the uncompressed pixels into their correct locations 
+        if (!zdata) return (FALSE);
+        if (!gfits_byteswap_zdata ((char *)zdata, Nzdata, 8)) return (FALSE);
+        if (!gfits_uncompress_data ((char *)zdata, Nzdata, cmptype, optname, optvalue, Noptions, out, &Nout, bytes_per_pixel)) return (FALSE);
+        if (!gfits_distribute_data (matrix, 8, out, Nout, otile, ztile, zscale, zzero)) return (FALSE);
+        // copy the uncompressed pixels into their correct locations        
         break;
       }
@@ -180 +190 @@
         short *zdata;
         zdata = gfits_varlength_column_pointer (ftable, &zdef, i, &Nzdata);
-        gfits_byteswap_zdata ((char *) zdata, Nzdata, 16);
-        gfits_uncompress_data ((char *) zdata, Nzdata, zdef.cmptype, optname, optvalue, Noptions, out, &Nout);
-        gfits_distribute_data (matrix, 16, out, Nout, otile, ztile, zscale, zzero);     // copy the uncompressed pixels into their correct locations 
+        if (!zdata) return (FALSE);
+        if (!gfits_byteswap_zdata ((char *)zdata, Nzdata, 16)) return (FALSE);
+        if (!gfits_uncompress_data ((char *)zdata, Nzdata, cmptype, optname, optvalue, Noptions, out, &Nout, bytes_per_pixel)) return (FALSE);
+        if (!gfits_distribute_data (matrix, 16, out, Nout, otile, ztile, zscale, zzero)) return (FALSE);
+        // copy the uncompressed pixels into their correct locations        
         break;
       }
@@ -188 +200 @@
       case 'J': {
         int *zdata;
-        zdata = gfits_varlength_column_pointer (ftable, &zdef, i, &Nzdata);
-        gfits_byteswap_zdata ((char *) zdata, Nzdata, 32);
-        gfits_uncompress_data ((char *) zdata, Nzdata, zdef.cmptype, optname, optvalue, Noptions, out, &Nout);
-        gfits_distribute_data (matrix, 32, out, Nout, otile, ztile, zscale, zzero);     // copy the uncompressed pixels into their correct locations 
+        if (!zdata) return (FALSE);
+        if (!gfits_byteswap_zdata ((char *)zdata, Nzdata, 32)) return (FALSE);
+        if (!gfits_uncompress_data ((char *)zdata, Nzdata, cmptype, optname, optvalue, Noptions, out, &Nout, bytes_per_pixel)) return (FALSE);
+        if (!gfits_distribute_data (matrix, 32, out, Nout, otile, ztile, zscale, zzero)) return (FALSE);
+        // copy the uncompressed pixels into their correct locations        
         break;
       }
 
       default:
-        abort();
+        fprintf (stderr, "invalid size for compressed data: %c\n", zdef.format);  
+        return (FALSE);
     }
 
@@ -212 +226 @@
 }
 
-int gfits_uncompress_data (char *zdata, int Nzdata, char *cmptype, char **optname, char **optvalue, int Nopt, char *outdata, int *Nout) {
-
-  // XXX allow lowercase value?
-  if (!strcasecmp(cmptype, "GZIP")) {
-    unsigned long tNout = *Nout;
-    uncompress ((Bytef *) outdata, &tNout, (Bytef *) zdata, Nzdata);
-    *Nout = tNout;
-    return (TRUE);
-  }
-
-  if (!strcasecmp(cmptype, "RICE")) {
+int gfits_tile_size (Matrix *matrix, int *otile, int *ztile) {
+
+  int i, Npixels, Ndimen;
+
+  // true sizes of this tile (in pixels)
+  Npixels = 1;
+  for (i = 0; i < matrix->Naxes; i++) {
+    Ndimen = MIN ((matrix->Naxis[i] - otile[i]*ztile[i]), ztile[i]);
+    Npixels *= Ndimen;
   }
   
-  if (!strcasecmp(cmptype, "PLIO")) {
-  }
-
-  if (!strcasecmp(cmptype, "HCOMPRESS")) {
-  }
-  return (TRUE);
+  return (Npixels);
 }
 
@@ -399 +406 @@
   return (TRUE);
 }
+
+int gfits_extension_is_compressed (Header *header) {
+
+    int has_extname, has_zimage, zimage;
+    char extname[80];
+
+    has_extname = gfits_scan (header, "EXTNAME", "%s", 1, extname);
+    has_zimage  = gfits_scan (header, "ZIMAGE", "%t", 1, &zimage);
+
+    if (has_zimage && zimage) return (TRUE);
+    if (has_extname) {
+        if (!strcmp (extname, "COMPRESSED_IMAGE")) return (TRUE);
+    }
+
+    return (FALSE);
+}
+
+int gfits_compressed_is_primary (Header *header) {
+
+    int has_ztension, has_zimage;
+    int ztension, zimage;
+
+    has_zimage   = gfits_scan (header, "ZIMAGE",   "%t", 1, &zimage);
+    has_ztension = gfits_scan (header, "ZTENSION", "%t", 1, &ztension);
+
+    if (has_zimage && zimage) return (TRUE);
+
+    return (FALSE);
+}

branches/eam_branch_20071015/Ohana/src/libfits/table/F_table_varlength.c

@@ -4 +4 @@
 int gfits_varlength_column_define (FTable *ftable, VarLengthColumn *def, int column) {
 
-  int i, N, Nv, Nb;
+  int i, Nv, Nb;
   char *p1, *p2, *p3;
   char field[81];
@@ -39 +39 @@
   if (*p1 == 'P') { return (FALSE); }
 
-  def[0].nbytes = 0;
-  if (*p1 == 'X') { def[0].nbytes = 1;  def[0].format = *p1; }
-  if (*p1 == 'L') { def[0].nbytes = 1;  def[0].format = *p1; }
-  if (*p1 == 'A') { def[0].nbytes = 1;  def[0].format = *p1; }
-  if (*p1 == 'B') { def[0].nbytes = 1;  def[0].format = *p1; }
-  if (*p1 == 'I') { def[0].nbytes = 2;  def[0].format = *p1; }
-  if (*p1 == 'J') { def[0].nbytes = 4;  def[0].format = *p1; }
-  if (*p1 == 'E') { def[0].nbytes = 4;  def[0].format = *p1; }
-  if (*p1 == 'D') { def[0].nbytes = 8;  def[0].format = *p1; }
-  if (*p1 == 'C') { def[0].nbytes = 8;  def[0].format = *p1; }
-  if (*p1 == 'M') { def[0].nbytes = 16; def[0].format = *p1; }
-  if (!def[0].nbytes) { return (FALSE); }
+  def->nbytes = 0;
+  if (*p1 == 'X') { def->nbytes = 1;  def->format = *p1; }
+  if (*p1 == 'L') { def->nbytes = 1;  def->format = *p1; }
+  if (*p1 == 'A') { def->nbytes = 1;  def->format = *p1; }
+  if (*p1 == 'B') { def->nbytes = 1;  def->format = *p1; }
+  if (*p1 == 'I') { def->nbytes = 2;  def->format = *p1; }
+  if (*p1 == 'J') { def->nbytes = 4;  def->format = *p1; }
+  if (*p1 == 'E') { def->nbytes = 4;  def->format = *p1; }
+  if (*p1 == 'D') { def->nbytes = 8;  def->format = *p1; }
+  if (*p1 == 'C') { def->nbytes = 8;  def->format = *p1; }
+  if (*p1 == 'M') { def->nbytes = 16; def->format = *p1; }
+  if (!def->nbytes) { return (FALSE); }
   
   /* scan columns to find column offset */
-  def[0].Nstart = 0;
-  for (i = 1; i < N; i++) {
+  def->Nstart = 0;
+  for (i = 1; i < column; i++) {
     sprintf (field, "TFORM%d", i);
     gfits_scan (ftable->header, field, "%s", 1, format);
     gfits_bintable_format (format, tmpline, &Nv, &Nb);
-    def[0].Nstart += Nv*Nb;
+    def->Nstart += Nv*Nb;
   }
 
-  if (!gfits_scan (ftable->header, "THEAP", "%d", 1, &def[0].theap)) {
-    def[0].theap = ftable->header->Naxis[0]*ftable->header->Naxis[1];
+  if (!gfits_scan (ftable->header, "THEAP", "%d", 1, &def->heap_start)) {
+    def->heap_start = ftable->header->Naxis[0]*ftable->header->Naxis[1];
   }
 
@@ -70 +70 @@
 void *gfits_varlength_column_pointer (FTable *ftable, VarLengthColumn *column, int row, int *length) {
 
-  int offset, Ny, *ptr;
+  int offset, Nx, *ptr;
 
   // find the values for the specified row
   // the values in the main table for this row and varlength column:
-  Ny = ftable->header->Naxis[1];
-  ptr = (int *) &ftable->buffer[row*Ny + column->Nstart];
-  offset = ptr[0];
-  *length = ptr[1];
+  Nx = ftable->header->Naxis[0];
+
+# define SWAP_WORD \
+  tmp = pchar[0]; pchar[0] = pchar[3]; pchar[3] = tmp; \
+  tmp = pchar[1]; pchar[1] = pchar[2]; pchar[2] = tmp;
+
+# ifdef BYTE_SWAP
+  char *pchar, tmp;
+  pchar = &ftable->buffer[row*Nx + column->Nstart];
+  SWAP_WORD;
+  pchar = &ftable->buffer[row*Nx + column->Nstart + 4];
+  SWAP_WORD;
+# endif
+
+  ptr = (int *) &ftable->buffer[row*Nx + column->Nstart];
+  *length = ptr[0];
+  offset = ptr[1];
 
   switch (column->format) {
@@ -85 +98 @@
     case 'B': {
       char *result;
-      result = (char *) ftable->buffer + column->theap + offset;
+      result = (char *) ftable->buffer + column->heap_start + offset;
       return result;
     }
     case 'I': {
       short *result;
-      result = (short *) ftable->buffer + column->theap + offset;
+      result = (short *) ftable->buffer + column->heap_start + offset;
       return result;
     }
     case 'J': {
       int *result;
-      result = (int *) ftable->buffer + column->theap + offset;
+      result = (int *) ftable->buffer + column->heap_start + offset;
       return result;
     }
@@ -101 +114 @@
     case 'C': {
       float *result;
-      result = (float *) ftable->buffer + column->theap + offset;
+      result = (float *) ftable->buffer + column->heap_start + offset;
       return result;
     }
@@ -107 +120 @@
     case 'M': {
       double *result;
-      result = (double *) ftable->buffer + column->theap + offset;
+      result = (double *) ftable->buffer + column->heap_start + offset;
       return result;
     }