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 38361) +++ /branches/eam_branches/ohana.20150429/src/libfits/matrix/F_compress_M.c (revision 38362) @@ -3,4 +3,6 @@ # include # define VERBOSE_DUMP 0 + +int gfits_collect_gzp2 (Matrix *matrix, char *raw, int Nraw, int raw_bitpix, int *otile, int oblank, int *ztile, int zblank, float zscale, float zzero); // the user needs to specify the various compression options: @@ -73,4 +75,11 @@ // define compression-specific keywords, update header as needed. if (!gfits_modify (theader, "ZCMPTYPE", "%s", 1, zcmptype)) ESCAPE; + + // define compression-specific keywords, update header as needed. + if (!strcasecmp(zcmptype, "GZIP_1") || !strcasecmp(zcmptype, "GZIP_2")) { + if ((header[0].bitpix == -32) || (header[0].bitpix == -64)) { + if (!gfits_modify (theader, "ZQUANTIZ", "%s", 1, "NONE")) ESCAPE; + } + } if (!gfits_modify_alt (theader, "ZSIMPLE", "%t", 1, TRUE)) ESCAPE; @@ -155,4 +164,12 @@ int cmp_pixsize = gfits_compressed_data_pixsize (zcmptype, header[0].bitpix, NULL, NULL, 0); + // define compression-specific keywords, update header as needed. + if (!strcasecmp(zcmptype, "RICE_1")) { + if (!gfits_modify (theader, "ZNAME1", "%s", 1, "BLOCKSIZE")) ESCAPE; + if (!gfits_modify (theader, "ZVAL1", "%d", 1, 32)) ESCAPE; + if (!gfits_modify (theader, "ZNAME2", "%s", 1, "BYTEPIX")) ESCAPE; + if (!gfits_modify (theader, "ZVAL2", "%d", 1, raw_pixsize)) ESCAPE; + } + // size of a pixel in the image tile -- this is probably not needed // int tile_pixsize = abs(header[0].bitpix) / 8; @@ -176,5 +193,9 @@ // copy the raw pixels from their native matrix locations to the temporary output buffer // 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 (!strcasecmp (zcmptype, "GZIP_2")) { + if (!gfits_collect_gzp2 (matrix, raw, Nraw, raw_bitpix, otile, oblank, ztile, zblank, zscale, zzero)) ESCAPE; + } else { + if (!gfits_collect_data (matrix, raw, Nraw, raw_bitpix, otile, oblank, ztile, zblank, zscale, zzero)) ESCAPE; + } if (VERBOSE_DUMP && (i == 0)) { @@ -214,7 +235,8 @@ } - // all other compression modes require swapping after compression - // (compresssion & decompression operate on native ENDIAN) - if (strcasecmp(zcmptype, "GZIP_1")) { + // GZIP_1 is swapped before compression; GZIP_2 implies swapping. All other + // compression modes require swapping after compression (compresssion & decompression + // operate on native ENDIAN) + if (strcasecmp(zcmptype, "GZIP_1") && strcasecmp(zcmptype, "GZIP_2")) { // Nzdata is number of bytes if (!gfits_byteswap_zdata (zdata, Nzdata, cmp_pixsize)) ESCAPE; @@ -383,2 +405,104 @@ return (TRUE); } + +// GZIP_2 organizes data in bytes by significance: (A1 A2 A3 A4), (B1 B2 B3 B4) -> (A1 B1), (A2, B2), (A3, B3), (A4, B4) +// NOTE that this organization depends on the machine endian state. +// raw_pixsize is the bytes / pixel for the input matrix +// place the raw image bytes for the current tile into the tile buffer +// * 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_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); // counter for current row in tile to copy + ALLOCATE (Ztile, int, matrix->Naxes); // 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]; + } + + // confirm the tile size: + assert (Nraw == Ztile[0]*Nline); + + // set the starting point of the tile in the matrix buffer + // 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)); + } + *rawptr = *srcptr; + 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); +} 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 38361) +++ /branches/eam_branches/ohana.20150429/src/libfits/matrix/F_compress_data.c (revision 38362) @@ -50,5 +50,7 @@ } - if (!strcasecmp(cmptype, "GZIP_1")) { + // GZIP_1 uses inflate / deflate, GZIP_2 does as well, but operates on a buffer with + // a different byte organization (gfits_distribute_gzp2) + if (!strcasecmp(cmptype, "GZIP_1") || !strcasecmp(cmptype, "GZIP_2")) { unsigned long Nbytes = Nrawpix * rawpix_size; 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 38361) +++ /branches/eam_branches/ohana.20150429/src/libfits/matrix/F_uncompress_M.c (revision 38362) @@ -3,4 +3,6 @@ # include # 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,5 +323,9 @@ // 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,5 +365,5 @@ 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,5 +380,5 @@ // number of lines in the tile (in pixels) - Nline = 1; + int Nline = 1; for (i = 1; i < matrix->Naxes; i++) { Nline *= Ztile[i]; @@ -385,12 +391,12 @@ // 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,2 +472,99 @@ } +// 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); +} + Index: /branches/eam_branches/ohana.20150429/src/libfits/matrix/F_uncompress_data.c =================================================================== --- /branches/eam_branches/ohana.20150429/src/libfits/matrix/F_uncompress_data.c (revision 38361) +++ /branches/eam_branches/ohana.20150429/src/libfits/matrix/F_uncompress_data.c (revision 38362) @@ -39,5 +39,5 @@ } - if (!strcasecmp(cmptype, "GZIP_1")) { + if (!strcasecmp(cmptype, "GZIP_1") || !strcasecmp(cmptype, "GZIP_2")) { // unsigned long tNout = *Nout * out_pixsize; // input value of *Nout is number of BYTES Index: /branches/eam_branches/ohana.20150429/src/libfits/test/compress.sh =================================================================== --- /branches/eam_branches/ohana.20150429/src/libfits/test/compress.sh (revision 38361) +++ /branches/eam_branches/ohana.20150429/src/libfits/test/compress.sh (revision 38362) @@ -2,6 +2,10 @@ macro test_image_all - foreach mode NONE GZIP_1 RICE_1 + # foreach mode NONE GZIP_1 RICE_1 + # foreach mode GZIP_1 + foreach mode GZIP_2 + # foreach mode RICE_1 foreach bitpix 8 16 32 -32 -64 + # foreach bitpix 32 echo ===== bitpix = $bitpix ===== test_image_single $bitpix $mode @@ -30,28 +34,28 @@ if ($cmpmode == NONE) return - if ($bitpix == -64) - echo "*** bitpix $bitpix fails fpack / funpack ***" - else - exec fpack -g -S test.raw.fits > test.fpk.fits - echo "load fpack version" - rd c test.fpk.fits -x 0 - set d = x - c - stat d - end - - if ($cmpmode == RICE_1) + if ($cmpmode == GZIP_1) + exec fpack -q 0 -g -S test.raw.fits > test.fpk.fits + end + if ($cmpmode == GZIP_2) + exec fpack -q 0 -g2 -S test.raw.fits > test.fpk.fits + end + if ($cmpmode == RICE_1) + exec fpack -r -S test.raw.fits > test.fpk.fits + end + echo "load fpack version" + rd c test.fpk.fits -x 0 + set d = x - c + stat d + + if ($cmpmode == RICE_1_0) echo "*** funpack fails on RICE_1 ****" return end - if ($bitpix == -64) - echo "*** bitpix $bitpix fails fpack / funpack ***" - else - echo "run funpack on our version" - exec funpack -S test.cmp.fits > test.fun.fits - rd e test.fun.fits - set d = x - e - stat d - end + echo "run funpack on our version" + exec funpack -S test.cmp.fits > test.fun.fits + rd e test.fun.fits + set d = x - e + stat d end Index: /branches/eam_branches/ohana.20150429/src/libfits/test/imagecomp.c =================================================================== --- /branches/eam_branches/ohana.20150429/src/libfits/test/imagecomp.c (revision 38361) +++ /branches/eam_branches/ohana.20150429/src/libfits/test/imagecomp.c (revision 38362) @@ -6,6 +6,6 @@ // char *cmptype[] = {"NONE", "GZIP_1", "GZIP_2", "PLIO_1", "RICE_1", "RICE_ONE", "HCOMPRESS_1", NULL}; -char *cmptype[] = {"NONE", "GZIP_1", "RICE_1", "RICE_ONE", NULL}; -int bitpix[] = {8, 16, 32, -32, -64, 0}; +// char *cmptype[] = {"NONE", "GZIP_1", "RICE_1", "RICE_ONE", NULL}; +// int bitpix[] = {8, 16, 32, -32, -64, 0}; // char *cmptype[] = {"NONE", "GZIP_1", NULL}; @@ -13,6 +13,6 @@ // int bitpix[] = {8, 16, 32, 0}; -// char *cmptype[] = {"RICE_1", NULL}; -// int bitpix[] = {-32, -64, 0}; +char *cmptype[] = {"GZIP_2", NULL}; +int bitpix[] = {-32, 0}; int main (int argc, char **argv) { @@ -43,5 +43,5 @@ } -# define NX 100 +# define NX 10 # define NY 10