Index: /branches/eam_branch_20071015/Ohana/src/libfits/Makefile =================================================================== --- /branches/eam_branch_20071015/Ohana/src/libfits/Makefile (revision 15352) +++ /branches/eam_branch_20071015/Ohana/src/libfits/Makefile (revision 15353) @@ -15,5 +15,5 @@ # programs may add their own internal requirements here -FULL_CFLAGS = $(BASE_CFLAGS) -fPIC +FULL_CFLAGS = $(BASE_CFLAGS) -fPIC -Wall -Werror FULL_CPPFLAGS = $(BASE_CPPFLAGS) FULL_LDFLAGS = $(BASE_LDFLAGS) @@ -32,5 +32,5 @@ $(HEAD)/F_copy_H.$(ARCH).o $(HEAD)/F_delete.$(ARCH).o \ $(HEAD)/F_read_XH.$(ARCH).o $(HEAD)/F_init_H.$(ARCH).o \ -$(HEAD)/version.$(ARCH).o +$(HEAD)/F_convert_H.$(ARCH).o $(HEAD)/version.$(ARCH).o MATRIX_OBJ = \ @@ -43,5 +43,5 @@ $(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_matrix.$(ARCH).o $(MATR)/F_compress_M.$(ARCH).o TABLE_OBJ = \ @@ -51,5 +51,6 @@ $(TABL)/F_define_column.$(ARCH).o $(TABL)/F_table_format.$(ARCH).o \ $(TABL)/F_set_column.$(ARCH).o $(TABL)/F_get_column.$(ARCH).o \ -$(TABL)/F_table_row.$(ARCH).o $(TABL)/F_free_T.$(ARCH).o +$(TABL)/F_table_row.$(ARCH).o $(TABL)/F_free_T.$(ARCH).o \ +$(TABL)/F_table_varlength.$(ARCH).o OBJS = $(HEADER_OBJ) $(MATRIX_OBJ) $(TABLE_OBJ) Index: /branches/eam_branch_20071015/Ohana/src/libfits/header/F_convert_H.c =================================================================== --- /branches/eam_branch_20071015/Ohana/src/libfits/header/F_convert_H.c (revision 15353) +++ /branches/eam_branch_20071015/Ohana/src/libfits/header/F_convert_H.c (revision 15353) @@ -0,0 +1,52 @@ +# include +# include + +int gfits_primary_to_extended (Header *header, char *exttype, char *comment) { + + int Ns, No; + char line[81]; + + // XXX check for valid exttype, comment, output string + snprintf (line, 81, "%-8s= '%-18s' / %-s ", "XTENSION", exttype, comment); + Ns = strlen (line); + No = 80 - Ns; + strncpy (header->buffer, line, Ns); + memset (&header->buffer[Ns], ' ', No); + + return (TRUE); +} + +// don't require the current to have SIMPLE +int gfits_modify_extended (Header *header, char *exttype, char *comment) { + + int Ns, No; + char line[81]; + + // XXX check for valid exttype, comment, output string + snprintf (line, 81, "%-8s= '%-18s' / %-s ", "XTENSION", exttype, comment); + Ns = strlen (line); + No = 80 - Ns; + strncpy (header->buffer, line, Ns); + memset (&header->buffer[Ns], ' ', No); + + return (TRUE); +} + +int gfits_extended_to_primary (Header *header, int simple, char *comment) { + + int Ns, No; + char line[81]; + + // XXX check for valid exttype, comment, output string + if (simple) { + snprintf (line, 81, "%-8s= %-18s T / %-s ", "SIMPLE", " ", comment); + } else { + snprintf (line, 81, "%-8s= %-18s F / %-s ", "SIMPLE", " ", comment); + } + Ns = strlen (line); + No = 80 - Ns; + strncpy (header->buffer, line, Ns); + memset (&header->buffer[Ns], ' ', No); + + return (TRUE); +} Index: /branches/eam_branch_20071015/Ohana/src/libfits/include/gfitsio.h =================================================================== --- /branches/eam_branch_20071015/Ohana/src/libfits/include/gfitsio.h (revision 15352) +++ /branches/eam_branch_20071015/Ohana/src/libfits/include/gfitsio.h (revision 15353) @@ -1,3 +1,5 @@ /* FITS specific macros and structures */ + +# include # ifndef GFITSIO @@ -91,5 +93,8 @@ int maxlen; int nbytes; + int Nstart; + int theap; char format; + char cmptype[81]; } VarLengthColumn; @@ -124,5 +129,9 @@ int gfits_vscan PROTO((Header *header, char *field, char *mode, int N, va_list argp)); int gfits_write_header PROTO((char *filename, Header *header)); -int gfits_matrix_size PROTO((Header *header)); +int gfits_data_size PROTO((Header *header)); +int gfits_extended_to_primary PROTO((Header *header, int simple, char *comment)); +int gfits_primary_to_extended PROTO((Header *header, char *exttype, char *comment)); +int gfits_modify_extended PROTO((Header *header, char *exttype, char *comment)); + /******************************* Matrix functions *************/ @@ -146,4 +155,8 @@ void gfits_set_matrix_value PROTO((Matrix *matrix, int x, int y, double value)); int gfits_write_matrix PROTO((char *filename, Matrix *matrix)); +Matrix *gfits_uncompress_image PROTO((Header *header, 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_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)); /******************************* Table functions *************/ @@ -188,3 +201,6 @@ int gfits_write_Theader PROTO((char *filename, Header *header)); +int gfits_varlength_column_define PROTO((FTable *ftable, VarLengthColumn *def, int column)); +void *gfits_varlength_column_pointer PROTO((FTable *ftable, VarLengthColumn *column, int row, int *length)); + #endif /* FITSIO */ Index: /branches/eam_branch_20071015/Ohana/src/libfits/matrix/F_compress_M.c =================================================================== --- /branches/eam_branch_20071015/Ohana/src/libfits/matrix/F_compress_M.c (revision 15352) +++ /branches/eam_branch_20071015/Ohana/src/libfits/matrix/F_compress_M.c (revision 15353) @@ -1,31 +1,56 @@ # include # include +# include # define MOD_KEYWORD(NAME,ZNAME,TYPE,IN,OUT) { \ gfits_scan (header, ZNAME, TYPE, 1, IN); \ gfits_delete_field (header, ZNAME, 1); \ - gfits_print (header, NAME, TYPE, 1, OUT); + gfits_print (header, NAME, TYPE, 1, OUT); } # define MOD_KEYWORD_REQUIRED(NAME,ZNAME,TYPE,IN,OUT) { \ - if (!gfits_scan (header, ZNAME, TYPE, 1, IN)) return (NULL); \ + if (!gfits_scan (header, ZNAME, TYPE, 1, IN)) ESCAPE; \ gfits_delete_field (header, ZNAME, 1); \ - gfits_print (header, NAME, TYPE, 1, OUT); - -Matrix *gfits_uncompress_image (Header *header, FTable *ftable) { - - int status; - int zimage; + gfits_print (header, NAME, TYPE, 1, OUT); } + +# 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 (NULL); } + +int gfits_uncompress_image (Header *header, Matrix *matrix, FTable *ftable, int primary) { + + int i, j, status, zimage, zcol, Ncur, Nzrows, Nout, Nzdata; char cmptype[80]; - char azis[10]; + char zaxis[10], naxis[10], key[10], word[81], exttype[81], checksum[81], datasum[81]; + int Noptions, NOPTIONS; + VarLengthColumn zdef; + float zscale, zzero; + + int *ztile = NULL; + int *otile = NULL; + int *ntile = NULL; + char **optname = NULL; + char **optvalue = NULL; + char *out = NULL; // is ZIMAGE present? status = gfits_scan (ftable->header, "ZIMAGE", "%t", 1, &zimage); - if (!status || !zimage) return NULL; + if (!status || !zimage) ESCAPE; // copy original header to output header - gfits_copy_header (&ftable->header, header); + gfits_copy_header (ftable->header, header); // extract compression-specific keywords, update header as needed. - if (!gfits_scan (header, "ZCMPTYPE", "%s", 1, cmptype)) return NULL; + if (!gfits_scan (header, "ZCMPTYPE", "%s", 1, cmptype)) ESCAPE; gfits_delete_field (header, "ZCMPTYPE", 1); @@ -35,6 +60,6 @@ for (i = 0; i < header->Naxes; i++) { snprintf (zaxis, 10, "ZNAXIS%d", i + 1); - snprintf (axis, 10, "NAXIS%d", i + 1); - MOD_KEYWORD_REQUIRED (zaxis, axis, "%d", &header->Naxis[i], header->Naxis[i]); + snprintf (naxis, 10, "NAXIS%d", i + 1); + MOD_KEYWORD_REQUIRED (zaxis, naxis, "%d", &header->Naxis[i], header->Naxis[i]); } @@ -48,7 +73,7 @@ } else { for (i = 1; i < header->Naxes; i++) { - snprintf (axis, 10, "ZTILE%d", i + 1); - if (!gfits_scan (header, axis, "%d", 1, &ztile[i])) return NULL; - gfits_delete_field (header, axis, 1); + snprintf (key, 10, "ZTILE%d", i + 1); + if (!gfits_scan (header, key, "%d", 1, &ztile[i])) ESCAPE; + gfits_delete_field (header, key, 1); } } @@ -59,12 +84,12 @@ ALLOCATE (optvalue, char *, NOPTIONS); for (Noptions = 0; TRUE; Noptions++) { - snprintf (axis, 10, "ZNAME%d", Noptions + 1); - if (!gfits_scan (header, axis, "%s", 1, word)) break; - gfits_delete_field (header, axis, 1); + snprintf (key, 10, "ZNAME%d", Noptions + 1); + if (!gfits_scan (header, key, "%s", 1, word)) break; + gfits_delete_field (header, key, 1); optname[Noptions] = strcreate (word); - snprintf (axis, 10, "ZVAL%d", Noptions + 1); - if (!gfits_scan (header, axis, "%s", 1, word)) ESCAPE; - gfits_delete_field (header, axis, 1); + snprintf (key, 10, "ZVAL%d", Noptions + 1); + if (!gfits_scan (header, key, "%s", 1, word)) ESCAPE; + gfits_delete_field (header, key, 1); optvalue[Noptions] = strcreate (word); @@ -75,13 +100,14 @@ } } + // XXX is Noptions set to the correct value? - // XXX get ZMASKCMP + // XXX can primary be defined on basis of the header? if (primary) { // XXX SIMPLE & XTENSION need to be set on the first line - gfits_scan (header, "ZSIMPLE", "%t", 1, &header->simple) return NULL; + if (!gfits_scan (header, "ZSIMPLE", "%t", 1, &header->simple)) ESCAPE; gfits_delete_field (header, "ZSIMPLE", 1); - gfits_print (header, "SIMPLE", "%d", 1, header->simple); + gfits_extended_to_primary (header, header->simple, "Image data"); MOD_KEYWORD ("ZEXTEND", "EXTEND", "%t", &header->extend, header->extend); @@ -89,7 +115,7 @@ } else { // XXX XTENSION needs to be set on the first line - gfits_scan (header, "ZTENSION", "%s", 1, exttype) return NULL; + if (!gfits_scan (header, "ZTENSION", "%s", 1, exttype)) ESCAPE; gfits_delete_field (header, "ZTENSION", 1); - gfits_print (header, "SIMPLE", "%d", 1, header->simple); + gfits_modify_extended (header, exttype, "Image extension"); MOD_KEYWORD ("ZPCOUNT", "PCOUNT", "%d", &header->pcount, header->pcount); @@ -99,22 +125,276 @@ MOD_KEYWORD ("ZHECKSUM", "ZHECKSUM", "%s", checksum, checksum); MOD_KEYWORD ("ZDATASUM", "DATASUM", "%s", datasum, datasum); + + zscale = 1; + gfits_scan (header, "ZSCALE", "%lf", 1, &zscale); + + zzero = 0; + gfits_scan (header, "ZZERO", "%lf", 1, &zzero); // find the COMPRESSED_DATA column (format should be 1PB, 1PI, 1PJ) for (i = 1; TRUE; i++) { snprintf (key, 10, "TTYPE%d", i); - if (!gfits_scan (&ftable->header, key, "%s", 1, ttype)) return NULL; - if (!strcmp (ttype, "COMPRESSED_DATA")) break; - } - compcol = i; - - if (!gfits_varlength_column_define (ftable, compdef, compcol)) return NULL; + if (!gfits_scan (ftable->header, key, "%s", 1, word)) ESCAPE; + if (!strcmp (word, "COMPRESSED_DATA")) break; + } + zcol = i; + + if (!gfits_varlength_column_define (ftable, &zdef, zcol)) ESCAPE; // create the output image - gfits_create_matrix (matrix, header); + 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); + 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]); + } // uncompress the data - - // ZSCALE & ZZERO convert uncompressed values to output values - // data = (uncompressed * ZSCALE) + ZZERO - + Nzrows = ftable->header->Naxis[1]; + for (i = 0; i < Nzrows; i++) { + 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 + + 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 + break; + } + + case 'I': { + 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 + break; + } + + 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 + break; + } + + default: + abort(); + } + + // update the tile counters, carrying to the next dimension if needed + for (j = 0; j < matrix->Naxes; j++) { + otile[j] ++; + if (otile[j] == ntile[j]) { + otile[j] = 0; + } else { + break; + } + } + } + return (TRUE); +} + +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")) { + } + if (!strcasecmp(cmptype, "PLIO")) { + } + + if (!strcasecmp(cmptype, "HCOMPRESS")) { + } + return (TRUE); +} + +// bitpix is the input data size/type +int gfits_distribute_data (Matrix *matrix, int bitpix, char *data, int Ndata, int *otile, int *ztile, float zscale, float zzero) { + + int i, j, start, offset, coord, Nline; + int *counter = NULL; + int *Ztile = NULL; + + // counter for current row in tile to copy + ALLOCATE (counter, int, matrix->Naxes); + memset (counter, 0, matrix->Naxes*sizeof(int)); + + // true sizes of this tile (in pixels) + ALLOCATE (Ztile, int, matrix->Naxes); + for (i = 0; i < matrix->Naxes; i++) { + Ztile[i] = MIN ((matrix->Naxis[i] - otile[i]*ztile[i]), ztile[i]); + } + + // number of lines in the tile (in pixels) + Nline = 1; + for (i = 1; i < matrix->Naxes; i++) { + Nline *= Ztile[i]; + } + + // double check reported size (needs the Isize) + assert (Ndata == 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] + ...)); + start = otile[matrix->Naxes-1]*ztile[matrix->Naxes-1]; + for (i = matrix->Naxes - 2; i >= 0; i--) { + coord = otile[i]*ztile[i]; + start = start*matrix->Naxis[i] + coord; + } + + // pixel offset in output array relative to tile start + offset = 0; + + // 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 +# define SCALE_AND_DIST(TYPE, SIZE) { \ + TYPE *Optr = (TYPE *) &matrix->buffer[SIZE*(offset + start)]; \ + for (j = 0; j < Ztile[0]; j++, Iptr++, Optr++) { \ + *Optr = *Iptr * zscale + zzero; \ + } } + + // this macro sets up the outer switch and calls above macro with all output bitpix options +# define SETUP_INSIZE(TYPE, SIZE) { \ + TYPE *Iptr = (TYPE *) &data[i*SIZE*Ztile[0]]; \ + switch (matrix->bitpix) { \ + case 8: \ + SCALE_AND_DIST (char, 1); \ + break; \ + case 16: \ + SCALE_AND_DIST (short, 2); \ + break; \ + case 32: \ + SCALE_AND_DIST (int, 4); \ + break; \ + case -32: \ + SCALE_AND_DIST (float, 4); \ + break; \ + case -64: \ + SCALE_AND_DIST (double, 8); \ + break; \ + default: \ + abort(); \ + } } + + // loop over lines in the tile + for (i = 0; i < Nline; i++) { + switch (bitpix) { + case 8: + SETUP_INSIZE (char, 1); + break; + case 16: + SETUP_INSIZE (short, 2); + break; + case 32: + SETUP_INSIZE (int, 4); + break; + default: + abort(); + } + + // 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 = 1; j < matrix->Naxes; j++) { + offset = offset*matrix->Naxis[j] + counter[j]; + } + } + + free (counter); + free (Ztile); + return (TRUE); +} + +int gfits_byteswap_zdata (char *zdata, int Nzdata, int bitpix) { + +# ifdef BYTE_SWAP + + int i; + char tmp; + + switch (bitpix) { + case +8: + case -8: + break; + + case +16: + case -16: + for (i = 0; i < Nzdata; i+=2) { + tmp = zdata[i]; + zdata[i] = zdata[i+1]; + zdata[i+1] = tmp; + } + break; + + case +32: + case -32: + for (i = 0; i < Nzdata; i+=4) { + tmp = zdata[i+1]; + zdata[i+1] = zdata[i+2]; + zdata[i+2] = tmp; + tmp = zdata[i]; + zdata[i] = zdata[i+3]; + zdata[i+3] = tmp; + } + break; + + case +64: + case -64: + for (i = 0; i < Nzdata; i+=8) { + tmp = zdata[i+0]; + zdata[i+0] = zdata[i+7]; + zdata[i+7] = tmp; + tmp = zdata[i+1]; + zdata[i+1] = zdata[i+6]; + zdata[i+6] = tmp; + tmp = zdata[i+2]; + zdata[i+2] = zdata[i+5]; + zdata[i+5] = tmp; + tmp = zdata[i+3]; + zdata[i+3] = zdata[i+4]; + zdata[i+4] = tmp; + } + break; + } +# endif + return (TRUE); +} Index: /branches/eam_branch_20071015/Ohana/src/libfits/table/F_table_format.c =================================================================== --- /branches/eam_branch_20071015/Ohana/src/libfits/table/F_table_format.c (revision 15352) +++ /branches/eam_branch_20071015/Ohana/src/libfits/table/F_table_format.c (revision 15353) @@ -194,5 +194,5 @@ int i, j, n, Nx, Ny, Nfields; int off, Nchar, Nval, Nbytes, status; - char *line, format[64], field[16], type[16]; + char format[64], field[16], type[16]; double tzero, tscale; char *tmpChar; @@ -233,5 +233,5 @@ for (j = 0; j < Ny; j++) { for (n = 0; n < Nval; n++) { - tmpChar = (int *)&ftable[0].buffer[j*Nx + n*Nbytes + off]; + tmpChar = (char *)&ftable[0].buffer[j*Nx + n*Nbytes + off]; *tmpChar = *tmpChar * tscale + tzero; } @@ -241,5 +241,5 @@ for (j = 0; j < Ny; j++) { for (n = 0; n < Nval; n++) { - tmpShort = (int *)&ftable[0].buffer[j*Nx + n*Nbytes + off]; + tmpShort = (short *)&ftable[0].buffer[j*Nx + n*Nbytes + off]; *tmpShort = *tmpShort * tscale + tzero; } @@ -264,5 +264,5 @@ int i, j, n, Nx, Ny, Nfields; int off, Nchar, Nval, Nbytes, status; - char *line, format[64], field[16], type[16]; + char format[64], field[16], type[16]; double tzero, tscale; char *tmpChar; @@ -303,5 +303,5 @@ for (j = 0; j < Ny; j++) { for (n = 0; n < Nval; n++) { - tmpChar = (int *)&ftable[0].buffer[j*Nx + n*Nbytes + off]; + tmpChar = (char *)&ftable[0].buffer[j*Nx + n*Nbytes + off]; *tmpChar = (*tmpChar - tzero) / tscale; } @@ -311,5 +311,5 @@ for (j = 0; j < Ny; j++) { for (n = 0; n < Nval; n++) { - tmpShort = (int *)&ftable[0].buffer[j*Nx + n*Nbytes + off]; + tmpShort = (short *)&ftable[0].buffer[j*Nx + n*Nbytes + off]; *tmpShort = (*tmpShort - tzero) / tscale; } Index: /branches/eam_branch_20071015/Ohana/src/libfits/table/F_table_varlength.c =================================================================== --- /branches/eam_branch_20071015/Ohana/src/libfits/table/F_table_varlength.c (revision 15352) +++ /branches/eam_branch_20071015/Ohana/src/libfits/table/F_table_varlength.c (revision 15353) @@ -4,10 +4,13 @@ int gfits_varlength_column_define (FTable *ftable, VarLengthColumn *def, int column) { - char field[80]; - char format[80]; + int i, N, Nv, Nb; + char *p1, *p2, *p3; + char field[81]; + char format[81]; + char tmpline[81]; // grab the value of TFORMn for this column snprintf (field, 80, "TFORM%d", column); - if (!gfits_scan (&ftable[0].header, field, "%s", 1, format)) return (FALSE); + if (!gfits_scan (ftable->header, field, "%s", 1, format)) return (FALSE); // find and remove the max field length element @@ -49,13 +52,67 @@ if (!def[0].nbytes) { return (FALSE); } + /* scan columns to find column offset */ + def[0].Nstart = 0; + for (i = 1; i < N; 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; + } + + if (!gfits_scan (ftable->header, "THEAP", "%d", 1, &def[0].theap)) { + def[0].theap = ftable->header->Naxis[0]*ftable->header->Naxis[1]; + } + return TRUE; } -void *gfits_varlength_column_pointer (FTable *ftable, VarLengthColumn *column, int row) { +void *gfits_varlength_column_pointer (FTable *ftable, VarLengthColumn *column, int row, int *length) { + + int offset, Ny, *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]; - - + switch (column->format) { + case 'X': + case 'L': + case 'A': + case 'B': { + char *result; + result = (char *) ftable->buffer + column->theap + offset; + return result; + } + case 'I': { + short *result; + result = (short *) ftable->buffer + column->theap + offset; + return result; + } + case 'J': { + int *result; + result = (int *) ftable->buffer + column->theap + offset; + return result; + } + case 'E': + case 'C': { + float *result; + result = (float *) ftable->buffer + column->theap + offset; + return result; + } + case 'D': + case 'M': { + double *result; + result = (double *) ftable->buffer + column->theap + offset; + return result; + } + case 'P': + default: + abort(); + } + abort(); } Index: /branches/eam_branch_20071015/Ohana/src/opihi/cmd.data/rd.c =================================================================== --- /branches/eam_branch_20071015/Ohana/src/opihi/cmd.data/rd.c (revision 15352) +++ /branches/eam_branch_20071015/Ohana/src/opihi/cmd.data/rd.c (revision 15353) @@ -86,4 +86,5 @@ } } + /* FITS extension by name */ if (ccdsel) { @@ -92,6 +93,6 @@ CCDKeyword = get_variable ("CCDKEYWORD"); if (CCDKeyword == (char *) NULL) { - gprint (GP_ERR, "CCDKEYWORD variable is not set; "); - gprint (GP_ERR, "using EXTNAME as default\n"); + // gprint (GP_ERR, "CCDKEYWORD variable is not set; "); + // gprint (GP_ERR, "using EXTNAME as default\n"); CCDKeyword = strcreate ("EXTNAME"); } @@ -106,5 +107,17 @@ return (FALSE); } - if (!gfits_scan (&buf[0].header, CCDKeyword, "%s", 1, ID)) { + + // for compressed data tables, EXTNAME may be duplicated, with the first one containing the + // word 'COMPRESSED_IMAGE'. in this case, check the second EXTNAME, if CCDKeyword == EXTNAME + // this may have to be a more obscure test specifically for 'imcopy' data... + // need to check each header, since file may contain a mix + + Nword = 1; + IsCompressed = FALSE; + if (gfits_extension_is_compressed_data (&buf[0].header)) { + if (!strcmp (CCDKeyword, "EXTNAME")) Nword = 2; + IsCompressed = TRUE; + } + if (!gfits_scan (&buf[0].header, CCDKeyword, "%s", Nword, ID)) { gprint (GP_ERR, "%s not in header\n", CCDKeyword); DeleteBuffer (buf); @@ -123,11 +136,7 @@ /* fix up header, if needed */ if (extend || ccdsel) { - static char simple[] = "SIMPLE = T / Standard FITS"; - int Ns, No; - - Ns = strlen (simple); - No = 80 - Ns; - strncpy (buf[0].header.buffer, simple, Ns); - memset (&buf[0].header.buffer[Ns], ' ', No); + if (!IsCompressed) { + gfits_extended_to_primary (buf[0].header, TRUE, "Standard FITS"); + } } else { gfits_fread_header (f, &buf[0].header); @@ -135,6 +144,7 @@ /* for JustHead, we skip reading the data segment */ + // XXX for compressed data, we need to convert the header to the equivalent uncompressed version if (JustHead) { - /* should be in CreateMatrix / CreateBuffer function */ + // XXX what are we doing here exactly? buf[0].header.Naxes = 0; ALLOCATE (buf[0].matrix.buffer, char, 1); @@ -164,7 +174,18 @@ /* load matrix data */ - sprintf (region, "-1 -1 -1 -1 %d %d", (plane - 1), plane); - status = gfits_fread_matrix_segment (f, &buf[0].matrix, &buf[0].header, region); - fclose (f); + if (IsCompressed) { + FTable ftable; + Header theader; + ftable.header = &theader; + gfits_copy_header (&buf[0].header, ftable.header); + status = gfits_fread_table (f, &ftable); // XXX does this do more than read the bytes? + status = gfits_uncompress_image (&buf[0].header, &buf[0].matrix, ftable, ccdsel || extend); + gfits_free_table (&ftable); + // XXX this currently does not work for a cube (we get a cube back, not a specific plane) + } else { + sprintf (region, "-1 -1 -1 -1 %d %d", (plane - 1), plane); + status = gfits_fread_matrix_segment (f, &buf[0].matrix, &buf[0].header, region); + fclose (f); + } if (!status) {