Index: /branches/eam_branches/ohana.20150429/src/libfits/include/gfitsio.h =================================================================== --- /branches/eam_branches/ohana.20150429/src/libfits/include/gfitsio.h (revision 38354) +++ /branches/eam_branches/ohana.20150429/src/libfits/include/gfitsio.h (revision 38355) @@ -192,7 +192,7 @@ int gfits_uncompress_image PROTO((Header *header, Matrix *matrix, FTable *ftable)); 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 oblank, int *ztile, int zblank, float zscale, float zzero)); - -int gfits_compress_data PROTO((char *zdata, int *Nzdata, char *cmptype, char **optname, char **optvalue, int Nopt, char *rawdata, int Nraw, int raw_pixsize)); +int gfits_distribute_data PROTO((Matrix *matrix, char *data, int Ndata, int bitpix, int *otile, int oblank, int *ztile, int zblank, float zscale, float zzero)); + +int gfits_compress_data PROTO((char *zdata, int *Nzdata, char *cmptype, char **optname, char **optvalue, int Nopt, char *rawdata, int Nraw, int raw_pixsize, int Nx, int Ny)); int gfits_compress_image PROTO((Header *header, Matrix *matrix, FTable *ftable, int *Ztile, char *zcmptype)); int gfits_collect_data PROTO((Matrix *matrix, char *raw, int Nraw, int raw_pixsize, int *otile, int oblank, int *ztile, int zblank, float zscale, float zzero)); @@ -205,9 +205,11 @@ int gfits_imtile_next PROTO((Matrix *matrix, int *ztile, int *ntile, int *otile)); -int gfits_byteswap_zdata PROTO((char *zdata, int Nzdata, int bitpix)); int gfits_extension_is_compressed_image PROTO((Header *header)); int gfits_extension_is_compressed_table PROTO((Header *header)); +int gfits_byteswap_zdata PROTO((char *zdata, int Nzdata, int pixsize)); +int gfits_compressed_data_pixsize PROTO((char *cmptype, int out_bitpix, char **optname, char **optvalue, int Noptions)); int gfits_uncompressed_data_pixsize PROTO((char *cmptype, int out_bitpix, char **optname, char **optvalue, int Noptions)); -int gfits_vartable_heap_pixsize PROTO((char format)); +int gfits_uncompressed_data_bitpix PROTO((char *cmptype, int out_bitpix, char **optname, char **optvalue, int Noptions)); +int gfits_vartable_heap_pixsize PROTO((char format)); int gfits_varlength_column_add_data PROTO((FTable *table, char *data, off_t Ndata, int row, VarLengthColumn *column)); Index: /branches/eam_branches/ohana.20150429/src/libfits/matrix/F_compress_M.c =================================================================== --- /branches/eam_branches/ohana.20150429/src/libfits/matrix/F_compress_M.c (revision 38354) +++ /branches/eam_branches/ohana.20150429/src/libfits/matrix/F_compress_M.c (revision 38355) @@ -18,5 +18,5 @@ // we do not try to support compression of an image with an associated table (not valid) -# define ESCAPE(A) { fprintf (stderr, "error in %s @ line %d\n", __func__, __LINE__); goto escape; } +# define ESCAPE { fprintf (stderr, "error in %s @ line %d\n", __func__, __LINE__); goto escape; } int gfits_compress_image (Header *header, Matrix *matrix, FTable *ftable, int *Ztile, char *zcmptype) { @@ -46,50 +46,50 @@ // creates an empty (Naxes = 2, Naxis[i] = 0) table header with XTENSION = BINTABLE, EXTNAME = COMPRESSED_IMAGE - if (!gfits_create_table_header (theader, "BINTABLE", "COMPRESSED_IMAGE")) ESCAPE(A); + if (!gfits_create_table_header (theader, "BINTABLE", "COMPRESSED_IMAGE")) ESCAPE; // by using "P" format, we are limited to 32bit pointers (2GB heap) // XXX how is the "B" format used? - if (!gfits_define_bintable_column (theader, "1PB(0)", "COMPRESSED_DATA", "compressed image data", "none", 1.0, 0.0)) ESCAPE (A); - if (!gfits_delete (theader, "TUNIT1", 1)) ESCAPE(A); + if (!gfits_define_bintable_column (theader, "1PB(0)", "COMPRESSED_DATA", "compressed image data", "none", 1.0, 0.0)) ESCAPE; + if (!gfits_delete (theader, "TUNIT1", 1)) ESCAPE; // allocates the default data array (ftable->buffer) - if (!gfits_create_table (theader, ftable)) ESCAPE (A); + if (!gfits_create_table (theader, ftable)) ESCAPE; // copy original header to output header (XXX this needs to be finished) - if (!gfits_copy_keywords_compress (header, theader)) ESCAPE (A); + if (!gfits_copy_keywords_compress (header, theader)) ESCAPE; // this allocates the pointer data array for the full set of tiles (NOT the heap) char *tmpbuffer = NULL; ALLOCATE_ZERO (tmpbuffer, char, Ntile); - if (!gfits_set_bintable_column (theader, ftable, "COMPRESSED_DATA", tmpbuffer, Ntile)) ESCAPE (A); + if (!gfits_set_bintable_column (theader, ftable, "COMPRESSED_DATA", tmpbuffer, Ntile)) ESCAPE; free (tmpbuffer); // add ZIMAGE from output header - if (!gfits_modify_alt (theader, "ZIMAGE", "%t", 1, TRUE)) ESCAPE (A); + if (!gfits_modify_alt (theader, "ZIMAGE", "%t", 1, TRUE)) ESCAPE; // define compression-specific keywords, update header as needed. - if (!gfits_modify (theader, "ZCMPTYPE", "%s", 1, zcmptype)) ESCAPE(A); - - if (!gfits_modify_alt (theader, "ZSIMPLE", "%t", 1, TRUE)) ESCAPE (A); + if (!gfits_modify (theader, "ZCMPTYPE", "%s", 1, zcmptype)) ESCAPE; + + if (!gfits_modify_alt (theader, "ZSIMPLE", "%t", 1, TRUE)) ESCAPE; // gfits_modify (header, "ZTENSION", "%s", 1, exttype); // supply the Z* header values from the source image header - if (!gfits_modify (theader, "ZBITPIX", "%d", 1, header[0].bitpix)) ESCAPE (A); - - if (!gfits_modify (theader, "ZNAXIS", "%d", 1, header[0].Naxes)) ESCAPE (A); + if (!gfits_modify (theader, "ZBITPIX", "%d", 1, header[0].bitpix)) ESCAPE; + + if (!gfits_modify (theader, "ZNAXIS", "%d", 1, header[0].Naxes)) ESCAPE; char keyword[11]; for (i = 0; i < header[0].Naxes; i++) { snprintf (keyword, 10, "ZNAXIS%d", i + 1); - if (!gfits_modify (theader, keyword, OFF_T_FMT, 1, header[0].Naxis[i])) ESCAPE (A); + if (!gfits_modify (theader, keyword, OFF_T_FMT, 1, header[0].Naxis[i])) ESCAPE; } for (i = 0; i < header->Naxes; i++) { snprintf (keyword, 10, "ZTILE%d", i + 1); - if (!gfits_modify (theader, keyword, "%d", 1, ztile[i])) ESCAPE(A); - } - - if (!gfits_modify (theader, "BSCALE", "%lf", 1, header[0].bscale)) ESCAPE (A); - if (!gfits_modify (theader, "BZERO", "%lf", 1, header[0].bzero)) ESCAPE (A); + if (!gfits_modify (theader, keyword, "%d", 1, ztile[i])) ESCAPE; + } + + if (!gfits_modify (theader, "BSCALE", "%lf", 1, header[0].bscale)) ESCAPE; + if (!gfits_modify (theader, "BZERO", "%lf", 1, header[0].bzero)) ESCAPE; // other keywords to define: @@ -128,15 +128,35 @@ VarLengthColumn zdef; - if (!gfits_varlength_column_define (ftable, &zdef, 1)) ESCAPE(A); - - // allocate the working buffers to the max tile size + if (!gfits_varlength_column_define (ftable, &zdef, 1)) ESCAPE; + + // pixel sizes and tile sizes: + // * We have Ntile tiles, each of Nx * Ny * Nz ... tile pixels + // -> A pixel in the tile (in the image) has size tile_pixsize = f(BITPIX) + // * We copy the pixels in the tile to a continuous buffer to be compressed. + // During this copy, we can scale the data (eg, floats -> ints), in which + // case the raw buffer may have a different pixel size (and type) than the + // image tile + // -> A pixel in the raw buffer is f(zcmptype,BITPIX): + // -- GZIP_1 : raw_pixsize = tile_pixsize (no scaling is performed) + // -- GZIP_2 : raw_pixsize = tile_pixsize (no scaling is performed) + // -- RICE_1 : raw pixels must be integers. raw_pixsize is user specified? (ZNAMEn,ZVALn : BYTEPIX = 1,2,4,8) + // -- PLIO_1 : raw_pixsize = 2 bytes + // -- HCOMPRESS_1 : raw pixels must be integers. + + // size (in pixels) of the largest tile int max_tile_size = gfits_imtile_maxsize (matrix, ztile); - // size of a pixel in the raw image tile: - int raw_pixsize = abs(header[0].bitpix) / 8; - - // size of a pixel in the raw image tile: - int tile_pixsize = gfits_uncompressed_data_pixsize (zcmptype, header[0].bitpix, NULL, NULL, 0); - + // size of a pixel in the image tile -- this is probably not needed + int tile_pixsize = abs(header[0].bitpix) / 8; + + // size of a pixel in the raw pixel buffer (before compression) + int raw_pixsize = gfits_uncompressed_data_pixsize (zcmptype, header[0].bitpix, NULL, NULL, 0); + int raw_bitpix = gfits_uncompressed_data_bitpix (zcmptype, header[0].bitpix, NULL, NULL, 0); + int cmp_pixsize = gfits_compressed_data_pixsize (zcmptype, header[0].bitpix, NULL, NULL, 0); + + fprintf (stderr, "raw_pixsize: %d, cmp_pixsize: %d, tile_pixsize: %d, raw_bitpix: %d\n", + raw_pixsize, cmp_pixsize, tile_pixsize, raw_bitpix); + + // allocate the buffer for compression work int Nzdata_alloc = raw_pixsize*max_tile_size + 100; ALLOCATE (raw, char, raw_pixsize*max_tile_size); @@ -144,5 +164,5 @@ // init the otile[] counters - if (!gfits_imtile_start (matrix, otile)) ESCAPE (A); + if (!gfits_imtile_start (matrix, otile)) ESCAPE; // compress the data : copy into a tile, compress the tile, then add to the output table @@ -150,25 +170,50 @@ for (i = 0; i < Ntile; i++) { - // size of the current tile + // size of the current tile in pixels int Nraw = gfits_tile_size (matrix, otile, ztile); // copy the raw pixels from their native matrix locations to the temporary output buffer - if (!gfits_collect_data (matrix, raw, Nraw, raw_pixsize, otile, oblank, ztile, zblank, zscale, zzero)) ESCAPE(A); - + // 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 (0 && (i == 0)) { + int k; + fprintf (stderr, "cmp mat: "); + for (k = 0; k < 32; k++) { fprintf (stderr, "0x%02hhx ", matrix->buffer[k]); } + fprintf (stderr, "\n"); + fprintf (stderr, "cmp raw: "); + for (k = 0; k < 32; k++) { fprintf (stderr, "0x%02hhx ", raw[k]); } + fprintf (stderr, "\n"); + } + + // gzip compresses bytes which are in BIG-ENDIAN order if (!strcasecmp(zcmptype, "GZIP_1")) { - if (!gfits_byteswap_zdata (raw, Nraw, raw_pixsize)) ESCAPE(A); - } - - // XXX the tile must not be > 2GB - - // optname, optvalue = NULL, Noptions = 0 + if (!gfits_byteswap_zdata (raw, Nraw, raw_pixsize)) ESCAPE; + } + if (0 && (i == 0)) { + int k; + fprintf (stderr, "cmp swp: "); + for (k = 0; k < 32; k++) { fprintf (stderr, "0x%02hhx ", raw[k]); } + fprintf (stderr, "\n"); + } + + // optname, optvalue = NULL, Noptions = 0 : defaults for RICE_1, HCOMPRESS_1 int Nzdata = Nzdata_alloc; // available space, replaced with actual output size on compression - if (!gfits_compress_data (zdata, &Nzdata, zcmptype, NULL, NULL, 0, raw, Nraw, tile_pixsize)) ESCAPE(A); - - if (!gfits_varlength_column_add_data (ftable, zdata, Nzdata, i, &zdef)) ESCAPE(A); - - // XXX need to activate - // if (FALSE) gfits_byteswap_zdata (zdata, Nzdata, zdata_pixsize) ESCAPE(A); + if (!gfits_compress_data (zdata, &Nzdata, zcmptype, NULL, NULL, 0, raw, Nraw, raw_pixsize, ztile[1], ztile[0])) ESCAPE; + if (0 && (i == 0)) { + int k; + fprintf (stderr, "cmp dat: "); + for (k = 0; k < 32; k++) { fprintf (stderr, "0x%02hhx ", zdata[k]); } + fprintf (stderr, "\n"); + } + + if (!gfits_varlength_column_add_data (ftable, zdata, Nzdata, i, &zdef)) ESCAPE; + + // all other compression modes require swapping after compression + // (compresssion & decompression operate on native ENDIAN) + if (strcasecmp(zcmptype, "GZIP_1")) { + if (!gfits_byteswap_zdata (zdata, Nzdata, cmp_pixsize)) ESCAPE; + } // update the otile[] counters, carrying to the next dimension if needed @@ -176,5 +221,11 @@ } - if (!gfits_varlength_column_finish (ftable, &zdef)) ESCAPE(A); + if (!gfits_varlength_column_finish (ftable, &zdef)) ESCAPE; + if (0) { + int k; + fprintf (stderr, "cmp tbl: "); + for (k = 0; k < 32; k++) { fprintf (stderr, "0x%02hhx ", ftable->buffer[k]); } + fprintf (stderr, "\n"); + } free (raw); @@ -209,5 +260,8 @@ // raw_pixsize is the bytes / pixel for the input matrix // place the raw image bytes for the current tile into the tile buffer -int gfits_collect_data (Matrix *matrix, char *raw, int Nraw, int raw_pixsize, int *otile, int oblank, int *ztile, int zblank, float zscale, float zzero) { +// * 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_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; @@ -246,66 +300,46 @@ // 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 - // note we copy Optr to Iptr -# define SCALE_AND_DIST(TYPE, SIZE) { \ - TYPE *Optr = (TYPE *) &matrix->buffer[SIZE*(offset + start)]; \ - for (j = 0; j < Ztile[0]; j++, Iptr++, Optr++) { \ - if (*Optr == oblank) { \ - *Iptr = zblank; \ +# define SCALE_AND_DIST_INT(TYPE, SIZE) { \ + TYPE *TILEptr = (TYPE *) &matrix->buffer[SIZE*(offset + start)]; \ + for (j = 0; j < Ztile[0]; j++, TILEptr++, RAWptr++) { \ + if (*TILEptr == oblank) { \ + *RAWptr = zblank; \ } else { \ - *Iptr = (*Optr - zzero) / zscale; \ + *RAWptr = (*TILEptr - zzero) / zscale; \ } } } // 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 - // note we copy Optr to Iptr -# define SCALE_AND_DIST_FLOAT(TYPE, SIZE) { \ - TYPE *Optr = (TYPE *) &matrix->buffer[SIZE*(offset + start)]; \ - for (j = 0; j < Ztile[0]; j++, Iptr++, Optr++) { \ - if (!isfinite(*Optr)) { \ - *Iptr = zblank; \ +# define SCALE_AND_DIST_FLOAT(TYPE, SIZE) { \ + TYPE *TILEptr = (TYPE *) &matrix->buffer[SIZE*(offset + start)]; \ + for (j = 0; j < Ztile[0]; j++, TILEptr++, RAWptr++) { \ + if (!isfinite(*TILEptr)) { \ + *RAWptr = zblank; \ } else { \ - *Iptr = (*Optr - zzero) / zscale; \ + *RAWptr = (*TILEptr - zzero) / zscale; \ } } } - // this macro sets up the outer switch and calls above macro with all output bitpix options -# define SETUP_INSIZE(TYPE, SIZE) { \ - TYPE *Iptr = (TYPE *) &raw[i*SIZE*Ztile[0]]; \ + // this macro sets up the outer switch and calls above the macro with all RAW buffer bitpix options + // +# define SETUP_RAWSIZE(TYPE, SIZE) { \ + TYPE *RAWptr = (TYPE *) &raw[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 (float, 4); \ - break; \ - case -64: \ - SCALE_AND_DIST_FLOAT (double, 8); \ - break; \ - default: \ - abort(); \ + case 8: SCALE_AND_DIST_INT (char, 1); break; \ + case 16: SCALE_AND_DIST_INT (short, 2); break; \ + case 32: SCALE_AND_DIST_INT (int, 4); break; \ + case -32: SCALE_AND_DIST_FLOAT (float, 4); break; \ + case -64: SCALE_AND_DIST_FLOAT (double, 8); break; \ + default: abort(); \ } } // loop over lines in the tile for (i = 0; i < Nline; i++) { - switch (raw_pixsize) { - case 1: - SETUP_INSIZE (char, 1); - break; - case 2: - SETUP_INSIZE (short, 2); - break; - case 4: - SETUP_INSIZE (int, 4); - break; - case 8: - SETUP_INSIZE (double, 8); - break; - default: - abort(); + switch (raw_bitpix) { + case 8: SETUP_RAWSIZE (char, 1); break; + case 16: SETUP_RAWSIZE (short, 2); break; + case 32: SETUP_RAWSIZE (int, 4); break; + case -32: SETUP_RAWSIZE (float, 4); break; + case -64: SETUP_RAWSIZE (double, 8); break; + default: abort(); } Index: /branches/eam_branches/ohana.20150429/src/libfits/matrix/F_compress_data.c =================================================================== --- /branches/eam_branches/ohana.20150429/src/libfits/matrix/F_compress_data.c (revision 38354) +++ /branches/eam_branches/ohana.20150429/src/libfits/matrix/F_compress_data.c (revision 38355) @@ -5,4 +5,6 @@ /* functions defined in ricecomp.c */ int fits_rcomp(int a[], int nx, unsigned char *c, int clen, int nblock); +int fits_rcomp_short(short a[], int nx, unsigned char *c, int clen, int nblock); +int fits_rcomp_byte(char a[], int nx, unsigned char *c, int clen, int nblock); int fits_rdecomp (unsigned char *c, int clen, unsigned int array[], int nx, int nblock); int fits_rdecomp_short (unsigned char *c, int clen, unsigned short array[], int nx, int nblock); @@ -29,8 +31,12 @@ # define ESCAPE(RET) { fprintf (stderr, "gzip error in %s @ %s:%d\n", __func__, __FILE__, __LINE__); return (RET); } +// XXX I'm putting Nx and Ny on the argument list -- only used by hcompress +// with a little work, these could go on the option list int gfits_compress_data (char *zdata, int *Nzdata, char *cmptype, char **optname, char **optvalue, int Nopt, - char *rawdata, int Nrawpix, int rawpix_size) { + char *rawdata, int Nrawpix, int rawpix_size, + int Nx, int Ny) { + int Nout; int status; @@ -59,5 +65,4 @@ } -# if (0) if (!strcasecmp(cmptype, "RICE_1")) { int i, blocksize; @@ -74,58 +79,57 @@ } - int Npix; - // Npix = *Nout * (out_pixsize / 4.0); - Npix = *Nout; - - Ninsum += Nzdata; - Noutsum += *Nout; + // Ninsum += Nzdata; + // Noutsum += *Nout; // fprintf (stderr, "%d comp bytes; %d uncomp 'pixels', totals: %d %d\n", Nzdata, Npix, Ninsum, Noutsum); - switch (out_pixsize) { + switch (rawpix_size) { case 4: - // rice decompression from the CFITSIO source tree : we need to tell it the expected number of pixels - // is also REQUIRES 4byte output, which is fairly stupid. - if (fits_rdecomp ((unsigned char *) zdata, Nzdata, (unsigned int *) outdata, Npix, blocksize)) { - fprintf (stderr, "error in rice decompression\n"); - return (FALSE); - } - return (TRUE); + // rice compression from the CFITSIO source tree + Nout = fits_rcomp ((int *) rawdata, Nrawpix, (unsigned char *) zdata, *Nzdata, blocksize); + break; case 2: - if (fits_rdecomp_short ((unsigned char *) zdata, Nzdata, (unsigned short *) outdata, Npix, blocksize)) { - fprintf (stderr, "error in rice decompression\n"); - return (FALSE); - } - return (TRUE); + Nout = fits_rcomp_short ((short *) rawdata, Nrawpix, (unsigned char *) zdata, *Nzdata, blocksize); + break; case 1: - if (fits_rdecomp_byte ((unsigned char *) zdata, Nzdata, (unsigned char *) outdata, Npix, blocksize)) { - fprintf (stderr, "error in rice decompression\n"); - return (FALSE); - } - return (TRUE); + Nout = fits_rcomp_byte ((char *) rawdata, Nrawpix, (unsigned char *) zdata, *Nzdata, blocksize); + break; default: - fprintf (stderr, "invalid output pixel size %d\n", out_pixsize); + fprintf (stderr, "invalid output pixel size %d\n", rawpix_size); + *Nzdata = 0; return (FALSE); } - + if (Nout < 0) { + fprintf (stderr, "error in rice decompression\n"); + *Nzdata = 0; + return (FALSE); + } + *Nzdata = Nout; + return TRUE; } if (!strcasecmp(cmptype, "PLIO_1")) { - int Npix; - Npix = pl_l2pi ((short *) zdata, 1, (int *) outdata, *Nout); - if (Npix != *Nout) { - fprintf (stderr, "error in plio decompression\n"); - return (FALSE); - } + // note the fortan starting point: zdata is decremented at start + Nout = pl_p2li ((int *) rawdata, 1, (short *) zdata, Nrawpix); + *Nzdata = Nout; + if (!Nout) return (FALSE); return (TRUE); } if (!strcasecmp(cmptype, "HCOMPRESS_1")) { - int Nx, Ny, scale; + long Nbytes = Nrawpix * rawpix_size; + int scale = 0; status = 0; // call hdecompress without smoothing - fits_hdecompress ((unsigned char *) zdata, FALSE, (int *) outdata, &Nx, &Ny, &scale, &status); + + if (rawpix_size == 4) { + fits_hcompress ((int *) rawdata, Nx, Ny, scale, (char *) zdata, &Nbytes, &status); + } + if (rawpix_size == 8) { + fits_hcompress64 ((long long *) rawdata, Nx, Ny, scale, (char *) zdata, &Nbytes, &status); + } + if (status) { fprintf (stderr, "error in hdecompress\n"); @@ -134,11 +138,7 @@ // fprintf (stderr, "decompression yields image %d x %d (scale: %d)\n", Nx, Ny, scale); - if (Nx * Ny != *Nout) { - fprintf (stderr, "error in hdecompress: mismatched output size\n"); - return (FALSE); - } + *Nzdata = Nbytes / rawpix_size; return (TRUE); } -# endif fprintf (stderr, "unknown compression %s\n", cmptype); Index: /branches/eam_branches/ohana.20150429/src/libfits/matrix/F_compress_utils.c =================================================================== --- /branches/eam_branches/ohana.20150429/src/libfits/matrix/F_compress_utils.c (revision 38354) +++ /branches/eam_branches/ohana.20150429/src/libfits/matrix/F_compress_utils.c (revision 38355) @@ -174,14 +174,19 @@ int i, Nbyte; - if (!strcasecmp(cmptype, "GZIP_1") || !strcasecmp(cmptype, "NONE")) { - if (out_bitpix == 8) return (1); - if (out_bitpix == 16) return (2); - if (out_bitpix == 32) return (4); + // GZIP_1, GZIP_2 use matched input & output pixels + if (!strcasecmp(cmptype, "GZIP_1") || + !strcasecmp(cmptype, "GZIP_2") || + !strcasecmp(cmptype, "NONE")) { + if (out_bitpix == 8) return (1); + if (out_bitpix == 16) return (2); + if (out_bitpix == 32) return (4); if (out_bitpix == -32) return (4); if (out_bitpix == -64) return (8); return (1); } - if (!strcasecmp(cmptype, "RICE_1")) { - + + // RICE_1 must operate on integer pixels + if (!strcasecmp(cmptype, "RICE_1") || + !strcasecmp(cmptype, "RICE_ONE")) { // if BYTEPIX option is specified, use that for Nbyte for (i = 0; i < Noptions; i++) { @@ -191,9 +196,15 @@ } } + if (out_bitpix == 8) return (1); + if (out_bitpix == 16) return (2); + if (out_bitpix == 32) return (4); return (4); } + + // PLIO_1 always results in 4-byte ints (not 2 byte ints?) if (!strcasecmp(cmptype, "PLIO_1")) { return (4); } + if (!strcasecmp(cmptype, "HCOMPRESS_1")) { if (out_bitpix == 8) return (4); @@ -204,4 +215,89 @@ } +int gfits_compressed_data_pixsize (char *cmptype, int out_bitpix, char **optname, char **optvalue, int Noptions) { + + int i, Nbyte; + + // GZIP_1, GZIP_2 should not be swapped after compression + if (!strcasecmp(cmptype, "GZIP_1") || + !strcasecmp(cmptype, "GZIP_2") || + !strcasecmp(cmptype, "NONE")) { + if (out_bitpix == 8) return (1); + if (out_bitpix == 16) return (1); + if (out_bitpix == 32) return (1); + if (out_bitpix == -32) return (1); + if (out_bitpix == -64) return (1); + return (1); + } + + // RICE_1 must operate on integer pixels + if (!strcasecmp(cmptype, "RICE_1") || + !strcasecmp(cmptype, "RICE_ONE")) { + // if BYTEPIX option is specified, use that for Nbyte + for (i = 0; i < Noptions; i++) { + if (!strcmp(optname[i], "BYTEPIX")) { + Nbyte = atoi (optvalue[i]); + return (Nbyte); + } + } + if (out_bitpix == 8) return (1); + if (out_bitpix == 16) return (2); + if (out_bitpix == 32) return (4); + return (4); + } + + // PLIO_1 always results in 4-byte ints (not 2 byte ints?) + if (!strcasecmp(cmptype, "PLIO_1")) { + return (4); + } + + if (!strcasecmp(cmptype, "HCOMPRESS_1")) { + if (out_bitpix == 8) return (4); + if (out_bitpix == 16) return (4); + return (8); + } + return (0); +} + +int gfits_uncompressed_data_bitpix (char *cmptype, int out_bitpix, char **optname, char **optvalue, int Noptions) { + + int i, Nbyte; + + // GZIP_1, GZIP_2 use matched input & output pixels + if (!strcasecmp(cmptype, "GZIP_1") || + !strcasecmp(cmptype, "GZIP_2") || + !strcasecmp(cmptype, "NONE")) { + return (out_bitpix); + } + + // RICE_1 must operate on integer pixels + if (!strcasecmp(cmptype, "RICE_1") || + !strcasecmp(cmptype, "RICE_ONE")) { + // if BYTEPIX option is specified, use that for Nbyte + for (i = 0; i < Noptions; i++) { + if (!strcmp(optname[i], "BYTEPIX")) { + Nbyte = atoi (optvalue[i]); + return (Nbyte*8); + } + } + if (out_bitpix == 8) return (out_bitpix); + if (out_bitpix == 16) return (out_bitpix); + if (out_bitpix == 32) return (out_bitpix); + return (32); + } + + // PLIO_1 always results in 4-byte ints (not 2 byte ints?) + if (!strcasecmp(cmptype, "PLIO_1")) { + return (32); + } + + if (!strcasecmp(cmptype, "HCOMPRESS_1")) { + if (out_bitpix == 8) return (32); + if (out_bitpix == 16) return (32); + return (64); + } + return (0); +} + int gfits_vartable_heap_pixsize (char format) { Index: /branches/eam_branches/ohana.20150429/src/libfits/matrix/F_uncompress_M.c =================================================================== --- /branches/eam_branches/ohana.20150429/src/libfits/matrix/F_uncompress_M.c (revision 38354) +++ /branches/eam_branches/ohana.20150429/src/libfits/matrix/F_uncompress_M.c (revision 38355) @@ -46,6 +46,4 @@ float zscale, zzero; - int zdata_pixsize, odata_pixsize; - int *ztile = NULL; int *otile = NULL; @@ -242,5 +240,5 @@ // size of an element in the vartable heap section - zdata_pixsize = gfits_vartable_heap_pixsize (zdef.format); + // zdata_pixsize = gfits_vartable_heap_pixsize (zdef.format); XXX // size of a pixel in the final image (not needed) @@ -248,10 +246,19 @@ // size of a pixel in the output from the decompression routine - odata_pixsize = gfits_uncompressed_data_pixsize (cmptype, header[0].bitpix, optname, optvalue, Noptions); - - int Nout_alloc = odata_pixsize*max_tile_size; + 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); + + int Nout_alloc = raw_pixsize*max_tile_size; ALLOCATE (out, char, Nout_alloc); off_t row; + + if (0) { + int k; + fprintf (stderr, "unc tbl: "); + for (k = 0; k < 32; k++) { fprintf (stderr, "0x%02hhx ", ftable->buffer[k]); } + fprintf (stderr, "\n"); + } // uncompress the data @@ -266,21 +273,46 @@ if (!zdata) ESCAPE; - // XXX not certain this is the correct place to do the swap (or if zdata_pixsize is - // the correct size to guide the swap) - if (strcasecmp(cmptype, "GZIP_1")) { - if (!gfits_byteswap_zdata (zdata, Nzdata, zdata_pixsize)) ESCAPE; + if (!gfits_byteswap_zdata (zdata, Nzdata, cmp_pixsize)) ESCAPE; + } + + if (0 && (row == 0)) { + int k; + fprintf (stderr, "unc dat: "); + for (k = 0; k < 32; k++) { fprintf (stderr, "0x%02hhx ", zdata[k]); } + fprintf (stderr, "\n"); } // gfits_uncompress_data uncompresses from zdata to the temporary output buffer which must be allocated // XXX the tile must not be > 2GB - if (!gfits_uncompress_data ((char *)zdata, Nzdata, cmptype, optname, optvalue, Noptions, out, &Nout, odata_pixsize)) ESCAPE; - + if (!gfits_uncompress_data ((char *)zdata, Nzdata, cmptype, optname, optvalue, Noptions, out, &Nout, raw_pixsize)) ESCAPE; + + if (0 && (row == 0)) { + int k; + fprintf (stderr, "unc swp: "); + for (k = 0; k < 32; k++) { fprintf (stderr, "0x%02hhx ", out[k]); } + fprintf (stderr, "\n"); + } + if (!strcasecmp(cmptype, "GZIP_1")) { - if (!gfits_byteswap_zdata (out, Nout, odata_pixsize)) ESCAPE; + if (!gfits_byteswap_zdata (out, Nout, raw_pixsize)) ESCAPE; + } + + if (0 && (row == 0)) { + int k; + fprintf (stderr, "unc raw: "); + for (k = 0; k < 32; k++) { fprintf (stderr, "0x%02hhx ", out[k]); } + fprintf (stderr, "\n"); } // copy the uncompressed pixels into their correct locations - if (!gfits_distribute_data (matrix, odata_pixsize, out, Nout, otile, oblank, ztile, zblank, zscale, zzero)) ESCAPE; + if (!gfits_distribute_data (matrix, out, Nout, raw_bitpix, otile, oblank, ztile, zblank, zscale, zzero)) ESCAPE; + + if (0 && (row == 0)) { + int k; + fprintf (stderr, "unc mat: "); + for (k = 0; k < 32; k++) { fprintf (stderr, "0x%02hhx ", matrix->buffer[k]); } + fprintf (stderr, "\n"); + } // update the tile counters, carrying to the next dimension if needed @@ -311,5 +343,5 @@ // bitpix is the input data size/type -int gfits_distribute_data (Matrix *matrix, int bitpix, char *data, int Ndata, int *otile, int oblank, int *ztile, int zblank, float zscale, float zzero) { +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; @@ -334,5 +366,5 @@ // double check reported size (needs the Isize) - assert (Ndata == Ztile[0]*Nline); + assert (Nraw == Ztile[0]*Nline); // set the starting point of the tile: @@ -350,53 +382,45 @@ // 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++) { \ - if (*Iptr == zblank) { \ - *Optr = oblank; \ - } else { \ - *Optr = *Iptr * zscale + zzero; \ +# define SCALE_AND_DIST_INT(TYPE, SIZE) { \ + TYPE *TILEptr = (TYPE *) &matrix->buffer[SIZE*(offset + start)]; \ + for (j = 0; j < Ztile[0]; j++, TILEptr++, RAWptr++) { \ + if (*TILEptr == zblank) { \ + *TILEptr = oblank; \ + } else { \ + *TILEptr = *RAWptr * zscale + zzero; \ } } } + // 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_FLOAT(TYPE, SIZE) { \ + TYPE *TILEptr = (TYPE *) &matrix->buffer[SIZE*(offset + start)]; \ + for (j = 0; j < Ztile[0]; j++, TILEptr++, RAWptr++) { \ + if (!isfinite(*TILEptr)) { \ + *TILEptr = oblank; \ + } else { \ + *TILEptr = *RAWptr * 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(); \ +# define SETUP_RAWSIZE(TYPE, SIZE) { \ + TYPE *RAWptr = (TYPE *) &raw[i*SIZE*Ztile[0]]; \ + switch (matrix->bitpix) { \ + case 8: SCALE_AND_DIST_INT (char, 1); break; \ + case 16: SCALE_AND_DIST_INT (short, 2); break; \ + case 32: SCALE_AND_DIST_INT (int, 4); break; \ + case -32: SCALE_AND_DIST_FLOAT (float, 4); break; \ + case -64: SCALE_AND_DIST_FLOAT (double, 8); break; \ + default: abort(); \ } } // loop over lines in the tile for (i = 0; i < Nline; i++) { - switch (bitpix) { - case 1: - SETUP_INSIZE (char, 1); - break; - case 2: - SETUP_INSIZE (short, 2); - break; - case 4: - SETUP_INSIZE (int, 4); - break; - case 8: - SETUP_INSIZE (double, 8); - break; - default: - abort(); + switch (raw_bitpix) { + case 8: SETUP_RAWSIZE (char, 1); break; + case 16: SETUP_RAWSIZE (short, 2); break; + case 32: SETUP_RAWSIZE (int, 4); break; + case -32: SETUP_RAWSIZE (float, 4); break; + case -64: SETUP_RAWSIZE (double, 8); break; + default: abort(); } Index: /branches/eam_branches/ohana.20150429/src/libfits/table/F_compress_T.c =================================================================== --- /branches/eam_branches/ohana.20150429/src/libfits/table/F_compress_T.c (revision 38354) +++ /branches/eam_branches/ohana.20150429/src/libfits/table/F_compress_T.c (revision 38355) @@ -161,5 +161,5 @@ int Nzdata = Nzdata_alloc; // available space, replaced with actual output size on compression - if (!gfits_compress_data (zdata, &Nzdata, fields[j].zctype, NULL, NULL, 0, raw, Nraw, fields[j].pixsize)) ESCAPE(A); + if (!gfits_compress_data (zdata, &Nzdata, fields[j].zctype, NULL, NULL, 0, raw, Nraw, fields[j].pixsize, 0, 0)) ESCAPE(A); if (!gfits_varlength_column_add_data (tgttable, zdata, Nzdata, i, &fields[j].zdef)) ESCAPE(A);