Index: trunk/Ohana/src/libfits/Makefile =================================================================== --- trunk/Ohana/src/libfits/Makefile (revision 38255) +++ trunk/Ohana/src/libfits/Makefile (revision 38441) @@ -1,5 +1,5 @@ default: install help: - @echo "make options: install libfits man clean dist" + @echo "make options: install libfits man clean dist test test.verbose test.clean" include ../../Makefile.System @@ -13,4 +13,6 @@ INC = $(HOME)/include MAN = $(HOME)/doc +TESTDIR = $(HOME)/test +TESTBIN = $(HOME)/test include ../../Makefile.Common @@ -19,4 +21,18 @@ FULL_CPPFLAGS = $(BASE_CPPFLAGS) -I$(EXT) FULL_LDFLAGS = $(BASE_LDFLAGS) -lohana + +TEST_CFLAGS = $(BASE_CFLAGS) +TEST_CPPFLAGS = $(BASE_CPPFLAGS) +TEST_LDFLAGS = $(BASE_LDFLAGS) -lFITS -lohana -ltap_ohana + + +TESTXTRA = tcomptiming +TESTPROG = zlib ricetest imagecomp tablecomp +$(TESTPROG) : % : $(TESTBIN)/% +$(TESTXTRA) : % : $(TESTBIN)/% +test: $(TESTPROG) + for i in $(TESTPROG); do $(TESTBIN)/$$i 2>&1 | $(TESTHARNESS); done +test.verbose: $(TESTPROG) + for i in $(TESTPROG); do $(TESTBIN)/$$i 2>&1 | $(TESTHARNESS) -v; done install: $(DESTLIB)/libFITS.a $(DESTLIB)/libFITS.$(DLLTYPE) $(DESTMAN)/fits.1 @@ -37,14 +53,16 @@ MATRIX_OBJ = \ -$(MATR)/F_add_M_value.$(ARCH).o $(MATR)/F_add_M.$(ARCH).o \ -$(MATR)/F_create_M.$(ARCH).o $(MATR)/F_divide_M.$(ARCH).o \ -$(MATR)/F_free_M.$(ARCH).o $(MATR)/F_copy_M.$(ARCH).o \ -$(MATR)/F_write_M.$(ARCH).o $(MATR)/F_get_M_value.$(ARCH).o \ -$(MATR)/F_multiply_M.$(ARCH).o $(MATR)/F_read_M.$(ARCH).o \ -$(MATR)/F_insert_M.$(ARCH).o $(MATR)/F_set_M_value.$(ARCH).o \ -$(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_uncompress_data.$(ARCH).o $(MATR)/F_swap_raw.$(ARCH).o +$(MATR)/F_add_M_value.$(ARCH).o $(MATR)/F_add_M.$(ARCH).o \ +$(MATR)/F_create_M.$(ARCH).o $(MATR)/F_divide_M.$(ARCH).o \ +$(MATR)/F_free_M.$(ARCH).o $(MATR)/F_copy_M.$(ARCH).o \ +$(MATR)/F_write_M.$(ARCH).o $(MATR)/F_get_M_value.$(ARCH).o \ +$(MATR)/F_multiply_M.$(ARCH).o $(MATR)/F_read_M.$(ARCH).o \ +$(MATR)/F_insert_M.$(ARCH).o $(MATR)/F_set_M_value.$(ARCH).o \ +$(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_utils.$(ARCH).o \ +$(MATR)/F_compress_M.$(ARCH).o $(MATR)/F_compress_data.$(ARCH).o \ +$(MATR)/F_uncompress_M.$(ARCH).o $(MATR)/F_uncompress_data.$(ARCH).o \ +$(MATR)/F_swap_raw.$(ARCH).o TABLE_OBJ = \ @@ -55,5 +73,6 @@ $(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_varlength.$(ARCH).o $(TABL)/F_copy_T.$(ARCH).o +$(TABL)/F_table_varlength.$(ARCH).o $(TABL)/F_copy_T.$(ARCH).o \ +$(TABL)/F_compress_T.$(ARCH).o $(TABL)/F_uncompress_T.$(ARCH).o EXTERN_OBJ = \ Index: trunk/Ohana/src/libfits/doc/notes-compress.txt =================================================================== --- trunk/Ohana/src/libfits/doc/notes-compress.txt (revision 38255) +++ trunk/Ohana/src/libfits/doc/notes-compress.txt (revision 38441) @@ -1,2 +1,20 @@ + +20150525 : + +uncompression steps: + + * read in a table: gfits_fread_ftable_data (f, &ftable, FALSE); (no byte swap) + * uncompress the image: gfits_uncompress_image (&buf[0].header, &buf[0].matrix, &ftable); + * get ZIMAGE from heaer & delete + * get ZCMPTYPE + * get ZBITPIX, ZNAXIS, ZNAXISnn + * get ZTILEn (if missing ZTILE1 = Nx, others 1) + * get ZNAME, ZVAL + * ZSIMPLE, ZTENSION, ZEXTEND, ZBLOCKED, ZPCOUNT, ZGCOUNT, etc + * ZSCALE, ZBLANK, OBLANK, ZZERO + + * data is stored in fits table COMPRESSED_DATA (1PB, 1PI, 1PJ) + +------------------------------------------------ TFORMn : rPt(e_max) Index: trunk/Ohana/src/libfits/extern/gzip.c =================================================================== --- trunk/Ohana/src/libfits/extern/gzip.c (revision 38255) +++ trunk/Ohana/src/libfits/extern/gzip.c (revision 38441) @@ -2,4 +2,5 @@ # include # include +# define EXTRA_VERBOSE 0 /** the two functions in this file re-implement the uncompress function of zlib. @@ -10,35 +11,35 @@ /* zlib.h -- interface of the 'zlib' general purpose compression library - version 1.2.3, July 18th, 2005 - - Copyright (C) 1995-2005 Jean-loup Gailly and Mark Adler - - This software is provided 'as-is', without any express or implied - warranty. In no event will the authors be held liable for any damages - arising from the use of this software. - - Permission is granted to anyone to use this software for any purpose, - including commercial applications, and to alter it and redistribute it - freely, subject to the following restrictions: - - 1. The origin of this software must not be misrepresented; you must not - claim that you wrote the original software. If you use this software - in a product, an acknowledgment in the product documentation would be - appreciated but is not required. - 2. Altered source versions must be plainly marked as such, and must not be - misrepresented as being the original software. - 3. This notice may not be removed or altered from any source distribution. - - Jean-loup Gailly Mark Adler - jloup@gzip.org madler@alumni.caltech.edu - - - The data format used by the zlib library is described by RFCs (Request for - Comments) 1950 to 1952 in the files http://www.ietf.org/rfc/rfc1950.txt - (zlib format), rfc1951.txt (deflate format) and rfc1952.txt (gzip format). + version 1.2.3, July 18th, 2005 + + Copyright (C) 1995-2005 Jean-loup Gailly and Mark Adler + + This software is provided 'as-is', without any express or implied + warranty. In no event will the authors be held liable for any damages + arising from the use of this software. + + Permission is granted to anyone to use this software for any purpose, + including commercial applications, and to alter it and redistribute it + freely, subject to the following restrictions: + + 1. The origin of this software must not be misrepresented; you must not + claim that you wrote the original software. If you use this software + in a product, an acknowledgment in the product documentation would be + appreciated but is not required. + 2. Altered source versions must be plainly marked as such, and must not be + misrepresented as being the original software. + 3. This notice may not be removed or altered from any source distribution. + + Jean-loup Gailly Mark Adler + jloup@gzip.org madler@alumni.caltech.edu + + + The data format used by the zlib library is described by RFCs (Request for + Comments) 1950 to 1952 in the files http://www.ietf.org/rfc/rfc1950.txt + (zlib format), rfc1951.txt (deflate format) and rfc1952.txt (gzip format). */ /* - The 'zlib' compression library provides in-memory compression and + The 'zlib' compression library provides in-memory compression and decompression functions, including integrity checks of the uncompressed data. This version of the library supports only one compression method @@ -46,5 +47,5 @@ stream interface. - Compression can be done in a single step if the buffers are large + Compression can be done in a single step if the buffers are large enough (for example if an input file is mmap'ed), or can be done by repeated calls of the compression function. In the latter case, the @@ -52,21 +53,21 @@ (providing more output space) before each call. - The compressed data format used by default by the in-memory functions is + The compressed data format used by default by the in-memory functions is the zlib format, which is a zlib wrapper documented in RFC 1950, wrapped around a deflate stream, which is itself documented in RFC 1951. - The library also supports reading and writing files in gzip (.gz) format + The library also supports reading and writing files in gzip (.gz) format with an interface similar to that of stdio using the functions that start with "gz". The gzip format is different from the zlib format. gzip is a gzip wrapper, documented in RFC 1952, wrapped around a deflate stream. - This library can optionally read and write gzip streams in memory as well. - - The zlib format was designed to be compact and fast for use in memory + This library can optionally read and write gzip streams in memory as well. + + The zlib format was designed to be compact and fast for use in memory and on communications channels. The gzip format was designed for single- file compression on file systems, has a larger header than zlib to maintain directory information, and uses a different, slower check method than zlib. - The library does not install any signal handler. The decoder checks + The library does not install any signal handler. The decoder checks the consistency of the compressed data, so the library should never crash even in case of corrupted input. @@ -116,37 +117,120 @@ } -// XXX ??? is cfitsio writing data which is incompatible with zlib ??? -int gfits_uncompress (Bytef *dest, uLongf *destLen, const Bytef *source, uLong sourceLen) -{ - z_stream stream; - int err; - - stream.next_in = (Bytef*)source; - stream.avail_in = (uInt)sourceLen; - /* Check for source > 64K on 16-bit machine: */ - if ((uLong)stream.avail_in != sourceLen) return Z_BUF_ERROR; - - stream.next_out = dest; - stream.avail_out = (uInt)*destLen; - if ((uLong)stream.avail_out != *destLen) return Z_BUF_ERROR; - - stream.zalloc = (alloc_func)0; - stream.zfree = (free_func)0; - - // MAX_WBITS = 15 - err = inflateInit2(&stream, -15); - if (err != Z_OK) return err; - - err = inflate(&stream, Z_FINISH); - if (err != Z_STREAM_END) { - inflateEnd(&stream); - if (err == Z_NEED_DICT || (err == Z_BUF_ERROR && stream.avail_in == 0)) - return Z_DATA_ERROR; - return err; - } - assert (stream.total_out <= *destLen); - *destLen = stream.total_out; - - err = inflateEnd(&stream); - return err; +# define ESCAPE(RET) { fprintf (stderr, "gzip error in %s @ %s:%d\n", __func__, __FILE__, __LINE__); return (RET); } + +int gfits_compress (Bytef *dest, uLongf *destLen, const Bytef *source, uLong sourceLen) { + + z_stream stream; + int i, err; + + stream.next_in = (Bytef*)source; + stream.avail_in = (uInt)sourceLen; + + /* Check for source > 64K on 16-bit machine: */ + if ((uLong)stream.avail_in != sourceLen) ESCAPE (Z_BUF_ERROR); + + stream.next_out = dest; + stream.avail_out = (uInt)*destLen; // allocated space + if ((uLong)stream.avail_out != *destLen) ESCAPE (Z_BUF_ERROR); + + stream.zalloc = Z_NULL; + stream.zfree = Z_NULL; + stream.opaque = Z_NULL; + + if (EXTRA_VERBOSE) { + fprintf (stderr, "inp cmp: "); + for (i = 0; i < sourceLen; i++) { + fprintf (stderr, "0x%02hhx ", source[i]); + } + fprintf (stderr, "\n"); + } + + // the '1' is the compression level: make this a user argument + // NOTE: cfitsio uses a fixed value of 1 + // I am using deflateInit2 because cfitsio expects gzip, not just zlib + // windowBits = 31 = (15 + 16) = (2^15 window bits) + (create a gzip stream) + err = deflateInit2(&stream, 1, Z_DEFLATED, 31, 8, Z_DEFAULT_STRATEGY); + if (0) fprintf (stderr, "inp buffers cmp 0: %d => %d => %d\n", stream.avail_in, stream.avail_out, (int) stream.total_out); + + if (err != Z_OK) ESCAPE(err); + + // XXX this is written to do the compression in a single pass. it could be re-done to have + // the compression occur in a series of steps + err = deflate(&stream, Z_FINISH); + if (err != Z_STREAM_END) { + err = deflateEnd(&stream); + ESCAPE (Z_BUF_ERROR); + } + + if (EXTRA_VERBOSE) { + fprintf (stderr, "out cmp: "); + for (i = 0; i < stream.total_out; i++) { + fprintf (stderr, "0x%02hhx ", dest[i]); + } + fprintf (stderr, "\n"); + } + + assert (stream.total_out <= *destLen); + *destLen = stream.total_out; + + err = deflateEnd(&stream); + return err; } + +// NOTE: CFITSIO uses gzip format, not just zlib +int gfits_uncompress (Bytef *dest, uLongf *destLen, const Bytef *source, uLong sourceLen) { + + z_stream stream; + int i, err; + + stream.next_in = (Bytef*)source; + stream.avail_in = (uInt)sourceLen; + + /* Check for source > 64K on 16-bit machine: */ + if ((uLong)stream.avail_in != sourceLen) ESCAPE (Z_BUF_ERROR); + + stream.next_out = dest; + stream.avail_out = (uInt)*destLen; + if ((uLong)stream.avail_out != *destLen) ESCAPE (Z_BUF_ERROR); + + stream.zalloc = Z_NULL; + stream.zfree = Z_NULL; + stream.opaque = Z_NULL; + + if (EXTRA_VERBOSE) { + fprintf (stderr, "inp unc: "); + for (i = 0; i < sourceLen; i++) { + fprintf (stderr, "0x%02hhx ", source[i]); + } + fprintf (stderr, "\n"); + } + + // windowBits = 47 = (15 + 32) = (2^15 window bits) + (test for either gzip or zlib + // streams). NOTE: modern versions of zlib (version > 1.2.3.4) allow for windowBits = 0 + // (or 16 or 32) to use the compressed stream windowBits value. IPP cluster still has + // zlib version = 1.2.3, so we cannot use this feature. + err = inflateInit2(&stream, 47); + if (err != Z_OK) ESCAPE (err); + + err = inflate(&stream, Z_FINISH); + if (err != Z_STREAM_END) { + inflateEnd(&stream); + if (err == Z_NEED_DICT || (err == Z_BUF_ERROR && stream.avail_in == 0)) ESCAPE (Z_DATA_ERROR); + ESCAPE(err); + } + + if (EXTRA_VERBOSE) { + fprintf (stderr, "out unc: "); + for (i = 0; i < stream.total_out; i++) { + fprintf (stderr, "0x%02hhx ", dest[i]); + } + fprintf (stderr, "\n"); + } + + assert (stream.total_out <= *destLen); + *destLen = stream.total_out; + + err = inflateEnd(&stream); + return err; +} + Index: trunk/Ohana/src/libfits/header/F_copy_H.c =================================================================== --- trunk/Ohana/src/libfits/header/F_copy_H.c (revision 38255) +++ trunk/Ohana/src/libfits/header/F_copy_H.c (revision 38441) @@ -30,2 +30,30 @@ return (TRUE); } + +int gfits_copy_header_ptr (Header *in, Header *out) { + + int i; + + if (!in) return FALSE; + if (!out) return FALSE; + + 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].Naxes = in[0].Naxes; + for (i = 0; i < FT_MAX_NAXES; i++) + out[0].Naxis[i] = in[0].Naxis[i]; + + out[0].datasize = in[0].datasize; + + out[0].pcount = in[0].pcount; + out[0].gcount = in[0].gcount; + out[0].bzero = in[0].bzero; + out[0].bscale = in[0].bscale; + + out[0].buffer = in[0].buffer; + + return (TRUE); +} Index: trunk/Ohana/src/libfits/header/F_modify.c =================================================================== --- trunk/Ohana/src/libfits/header/F_modify.c (revision 38255) +++ trunk/Ohana/src/libfits/header/F_modify.c (revision 38441) @@ -48,5 +48,5 @@ strncpy (comment, qe, p + 80 - qe); } - gfits_pad_ending (comment, ' ', 82); /* comment must contain spaces to the end */ + gfits_pad_ending (comment, 0x20, 82); /* comment must contain spaces to the end */ /* write the numeric modes */ @@ -134,5 +134,5 @@ strncpy (comment, qe, p + 80 - qe); } - gfits_pad_ending (comment, ' ', 82); /* comment must contain spaces to the end */ + gfits_pad_ending (comment, 0x20, 82); /* comment must contain spaces to the end */ /* write the boolean mode */ @@ -166,5 +166,5 @@ if (!ptr) goto invalid; strncpy (comment, ptr, 80); - gfits_pad_ending (comment, ' ', 82); /* comment must contain spaces to the end */ + gfits_pad_ending (comment, 0x20, 82); /* comment must contain spaces to the end */ snprintf (string, 81, "%-8s= %s / %-s", field, data, comment); /* this will keep the original line, but truncate the comment */ @@ -296,15 +296,23 @@ */ +// XXX this is an absurd patch on a gcc bug: memset has a problem is value != 0 and N is constant +void myMemset (char *ptr, int value, size_t N) { + + char *p = ptr; + + size_t i = 0; + for (i = 0; i < N; i++, p++) { + *p = value; + } +} + /* fill 'line' with Nbyte space from first NULL to last byte with value */ void gfits_pad_ending (char *line, char value, int Nbyte) { - - char *p; - int N; line[Nbyte-1] = 0; - p = line + strlen (line); - N = MAX (Nbyte - strlen (line) - 1, 0); - memset (p, value, N); + char *p = line + strlen (line); + size_t N = MAX (Nbyte - strlen (line) - 1, 0); + myMemset (p, value, N); } Index: trunk/Ohana/src/libfits/include/gfitsio.h =================================================================== --- trunk/Ohana/src/libfits/include/gfitsio.h (revision 38255) +++ trunk/Ohana/src/libfits/include/gfitsio.h (revision 38441) @@ -73,4 +73,5 @@ off_t datasize; // size of the buffer (including block padding at the end) off_t validsize; // size of the valid portion of the table (< validsize if file is short) + off_t heap_start; // byte offset to start of HEAP } FTable; @@ -85,12 +86,32 @@ } VTable; -typedef struct { - int maxlen; // max size of all table rows - int nbytes; // number of bytes per column element - int Nstart; // byte offset of this column relative to first column - off_t heap_start; // byte offset to start of HEAP - char format; // data format character (one of: XLABIJEDCM) - char mode; // access data pointer type (one of P or Q) +// description / metadata needed for a single varlength column. A varlength column has a +// 2-element array giving the (length,offset) of data in the heap. Below, 'real data' are +// the values in the heap; 'metadata' are the values in the main table. +typedef struct { + int column; // field number (of the metadata) + char format; // format character for real data (one of: XLABIJKEDCM) + char mode; // format character for metadata (one of P [float] or Q [double]) + int maxlen; // max size of all table rows for this column + int nbytes; // number of bytes per real data column element + int offset; // byte offset of the metadata column relative to first column } VarLengthColumn; + +typedef struct { + char ttype[16]; // TTYPE field of original table + char ttype_cmt[80]; // comment associated with TTYPE + char tunit[16]; // TUNIT field of original table + char tunit_cmt[80]; // comment associated with TUNIT + char tformat[16]; // TFORM field of original table + char tformat_cmt[80]; // comment associated with TFORM + char zctype[80]; // compression type for this field + VarLengthColumn zdef; // description of the output variable length column + + char datatype[80]; // named data type associated with ttype (eg, 18J -> int, 18D -> double) + int Nvalues; // number of values per row (eg 18J -> 18 int values) + int pixsize; // number of bytes per value of this row (eg, 18J -> 4 bytes) + int rowsize; // number of bytes for a full row (Nvalues * pixsize) + int offset; // byte offset of the real data column relative to first column +} TableField; # ifndef PROTO @@ -106,4 +127,5 @@ char *gfits_keyword_end PROTO((char *line)); int gfits_copy_header PROTO((Header *in, Header *out)); +int gfits_copy_header_ptr PROTO((Header *in, Header *out)); int gfits_create_header PROTO((Header *header)); int gfits_delete PROTO((Header *header, char *field, int N)); @@ -117,6 +139,6 @@ int gfits_modify PROTO((Header *header, char *field, char *mode, int N,...)) OHANA_FORMAT(printf, 3, 5); int gfits_print PROTO((Header *header, char *field, char *mode, int N,...)) OHANA_FORMAT(printf, 3, 5); -int gfits_modify_alt PROTO((Header *header, char *field, char *mode, int N,...)); -int gfits_print_alt PROTO((Header *header, char *field, char *mode, int N,...)); +int gfits_modify_alt PROTO((Header *header, char *field, char *mode, int N,...)); // do not use a FORMAT: non-standard fmt chars +int gfits_print_alt PROTO((Header *header, char *field, char *mode, int N,...)); // do not use a FORMAT: non-standard fmt chars int gfits_find_Xheader PROTO((FILE *f, Header *header, char *extname)); int gfits_read_Xheader PROTO((char *filename, Header *header, int N)); @@ -124,5 +146,5 @@ int gfits_save_header PROTO((FILE *f, Header *header)); int gfits_scan PROTO((Header *header, char *field, char *mode, int N,...)) OHANA_FORMAT(scanf, 3, 5); -int gfits_scan_alt PROTO((Header *header, char *field, char *mode, int N,...)); +int gfits_scan_alt PROTO((Header *header, char *field, char *mode, int N,...)); // do not use a FORMAT: non-standard fmt chars int gfits_set_unsign_mode PROTO((int mode)); int gfits_stripwhite PROTO((char *string)); @@ -133,4 +155,5 @@ off_t gfits_data_min_size PROTO((Header *header)); off_t gfits_data_pad_size PROTO((off_t rawsize)); +off_t gfits_heap_start 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)); @@ -166,13 +189,36 @@ void gfits_set_matrix_value PROTO((Matrix *matrix, off_t x, off_t y, double value)); int gfits_write_matrix PROTO((char *filename, Matrix *matrix)); + 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_byteswap_zdata PROTO((char *zdata, int Nzdata, int bitpix)); +int gfits_uncompress_data PROTO((char *zdata, int Nzdata, char *cmptype, char **optname, char **optvalue, int Nopt, char *outdata, int *Nout, int Nout_alloc, int out_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)); + +int gfits_copy_keywords_compress PROTO((Header *srchead, Header *tgthead)); int gfits_swap_raw PROTO((Matrix *matrix)); -int gfits_extension_is_compressed PROTO((Header *header)); + off_t gfits_tile_size PROTO((Matrix *matrix, int *otile, int *ztile)); +int gfits_imtile_maxsize PROTO((Matrix *matrix, int *ztile)); +int gfits_imtile_count PROTO((Matrix *matrix, int *ztile, int *ntile)); +int gfits_imtile_start PROTO((Matrix *matrix, int *otile)); +int gfits_imtile_next PROTO((Matrix *matrix, int *ztile, int *ntile, int *otile)); + +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_cmptype_valid PROTO((char *cmptype)); + +int gfits_varlength_column_add_data PROTO((FTable *table, char *data, off_t Ndata, int row, VarLengthColumn *column)); +void *gfits_varlength_column_pointer PROTO((FTable *ftable, VarLengthColumn *column, off_t row, off_t *length)); +int gfits_varlength_column_define PROTO((FTable *ftable, VarLengthColumn *def, int column)); +int gfits_varlength_column_finish PROTO((FTable *ftable, VarLengthColumn *def)); +int gfits_byteswap_varlength_column PROTO((FTable *ftable, int column)); /******************************* Table functions *************/ @@ -183,4 +229,5 @@ char *gfits_table_print PROTO((FTable *ftable,...)); int gfits_add_rows PROTO((FTable *ftable, char *data, off_t Nrow, off_t Nbytes)); +int gfits_set_table_rows PROTO((Header *header, FTable *table, off_t Nrows)); int gfits_bintable_format PROTO((char *format, char *type, int *Nval, int *Nbytes)); int gfits_create_table PROTO((Header *header, FTable *ftable)); @@ -199,7 +246,9 @@ int gfits_fwrite_ftable_range PROTO((FILE *f, FTable *table, off_t start, off_t Nrows, off_t Ndisk, off_t Ntotal)); int gfits_get_bintable_column PROTO((Header *header, FTable *table, char *label, void **data)); +int gfits_get_bintable_column_raw PROTO((Header *header, FTable *table, char *label, void **data, char nativeOrder)); int gfits_get_bintable_column_type PROTO((Header *header, char *label, char *type, int *Nval)); +void *gfits_get_bintable_column_data PROTO((Header *header, FTable *table, char *label, char *type, off_t *Nrow, int *Ncol)); int gfits_get_bintable_column_type_by_N PROTO((Header *header, int N, char *type, int *Nval)); -void *gfits_get_bintable_column_data PROTO((Header *header, FTable *table, char *label, char *type, off_t *Nrow, int *Ncol)); +void *gfits_get_bintable_column_data_raw PROTO((Header *header, FTable *table, char *label, char *type, off_t *Nrow, int *Ncol, char nativeOrder)); int gfits_get_table_column PROTO((Header *header, FTable *table, char *label, void **data)); int gfits_get_table_column_type PROTO((Header *header, char *label, char *type)); @@ -207,5 +256,5 @@ int gfits_read_table PROTO((char *filename, FTable *ftable)); int gfits_set_bintable_column PROTO((Header *header, FTable *table, char *label, void *data, off_t Nrow)); -int gfits_set_bintable_column_reformat PROTO((Header *header, FTable *table, char *label, char *intype, void *data, off_t Nrow)); +int gfits_set_bintable_column_reformat PROTO((Header *header, FTable *table, char *label, char *intype, void *data, off_t Nrow, char nativeOrder)); int gfits_set_table_column PROTO((Header *header, FTable *table, char *label, void *data, off_t Nrow)); int gfits_table_column PROTO((FTable *ftable, char *field, char *mode,...)) OHANA_FORMAT(printf, 3, 4); @@ -219,4 +268,5 @@ int gfits_copy_ftable PROTO((FTable *in, FTable *out)); int gfits_copy_vtable PROTO((VTable *in, VTable *out)); +int gfits_copy_ftable_ptr PROTO((FTable *in, FTable *out)); int gfits_create_Theader PROTO((Header *header, char *type)); @@ -227,6 +277,11 @@ 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, off_t row, off_t *length)); +int gfits_compress_table PROTO((FTable *srctable, FTable *tgttable, int ztilelen, char *zcmptype)); +int gfits_uncompress_table PROTO((FTable *srctable, FTable *tgttable)); + +// XXX EAM for testing +int gfits_dump_raw_table (FTable *table, char *message); +int gfits_dump_cmp_table (FTable *table, char *message); + #endif /* FITSIO */ Index: trunk/Ohana/src/libfits/matrix/F_compress_M.c =================================================================== --- trunk/Ohana/src/libfits/matrix/F_compress_M.c (revision 38255) +++ trunk/Ohana/src/libfits/matrix/F_compress_M.c (revision 38441) @@ -2,313 +2,306 @@ # include # include - -# 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 (FALSE); } - -# define MOD_KEYWORD(ZNAME,NAME,TYPE,IN,OUT) { \ - if (gfits_scan (header, ZNAME, TYPE, 1, IN)) { \ - gfits_modify (header, NAME, TYPE, 1, OUT); \ - } \ - gfits_delete (header, ZNAME, 1); } - -# define MOD_KEYWORD_ALT(ZNAME,NAME,TYPE,IN,OUT) { \ - if (gfits_scan_alt (header, ZNAME, TYPE, 1, IN)) { \ - gfits_modify_alt (header, NAME, TYPE, 1, OUT); \ - } \ - gfits_delete (header, ZNAME, 1); } - -# define MOD_KEYWORD_REQUIRED(ZNAME,NAME,TYPE,IN,OUT) { \ - if (!gfits_scan (header, ZNAME, TYPE, 1, IN)) ESCAPE; \ - gfits_delete (header, ZNAME, 1); \ - gfits_modify (header, NAME, TYPE, 1, OUT); } - -int gfits_uncompress_image (Header *header, Matrix *matrix, FTable *ftable) { - - off_t Nzdata, Nzrows, zcol; - int i, j, status, zimage, Nout, max_tile_size; - char cmptype[80]; - char zaxis[10], naxis[10], key[10], word[81], exttype[81], checksum[81], datasum[81]; - int Noptions, NOPTIONS, zblank, oblank; - VarLengthColumn zdef; - float zscale, zzero; - - int zdata_pixsize, odata_pixsize; - - int *ztile = NULL; +# 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: +// ztile[i] (may be NULL) +// zcmptype + +// steps to convert an image to a compressed table: +// * determine the number of tiles (from ztile[] and image size) +// * construct an empty table with the right dimensions (1 column, Ntile rows) +// * loop over tiles +// * extract tile data to a buffer +// * compress the buffer data +// * insert in table heap +// * update header keywords + +// we do not try to support compression of an image with an associated table (not valid) + +# 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) { + + int i; + + // pointers to track the tile data int *otile = NULL; int *ntile = NULL; - char **optname = NULL; - char **optvalue = NULL; - char *out = NULL; + char *raw = NULL; char *zdata = NULL; - Noptions = 0; - - // is ZIMAGE present? - // NOTE target of %t must be int length - status = gfits_scan_alt (ftable->header, "ZIMAGE", "%t", 1, &zimage); - if (!status || !zimage) ESCAPE; - - // copy original header to output header - gfits_copy_header (ftable->header, header); - - // delete ZIMAGE from output header - gfits_delete (header, "ZIMAGE", 1); - - // extract compression-specific keywords, update header as needed. - if (!gfits_scan (header, "ZCMPTYPE", "%s", 1, cmptype)) ESCAPE; - gfits_delete (header, "ZCMPTYPE", 1); - - MOD_KEYWORD_REQUIRED ("ZBITPIX", "BITPIX", "%d", &header->bitpix, header->bitpix); - MOD_KEYWORD_REQUIRED ("ZNAXIS", "NAXIS", "%d", &header->Naxes, header->Naxes); + int *ztile = Ztile; + + Header *theader = ftable->header; + + if (!gfits_cmptype_valid(zcmptype)) ESCAPE; + + // determine the number of tiles (from ztile[] and image size) + if (!ztile) { + ALLOCATE (ztile, int, matrix->Naxes); + ztile[0] = matrix->Naxis[0]; + for (i = 1; i < matrix->Naxes; i++) ztile[i] = 1; + } + ALLOCATE (ntile, int, matrix->Naxes); + ALLOCATE (otile, int, matrix->Naxes); + + int Ntile = gfits_imtile_count (matrix, ztile, ntile); + + // 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; + + // 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; + if (!gfits_delete (theader, "TUNIT1", 1)) ESCAPE; + + // allocates the default data array (ftable->buffer) + 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; + + // 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; + free (tmpbuffer); + + // add ZIMAGE from output header + 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; + + // 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; + // 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; + + 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; + } for (i = 0; i < header->Naxes; i++) { - snprintf (zaxis, 10, "ZNAXIS%d", i + 1); - snprintf (naxis, 10, "NAXIS%d", i + 1); - MOD_KEYWORD_REQUIRED (zaxis, naxis, OFF_T_FMT, &header->Naxis[i], header->Naxis[i]); - } - - // set up the tile sizes : ztile is the default size of the tile in the Nth dimension - // the actual tile size may be smaller at the edge of a dimension. if the ZTILEn - // entries are not found, default to [Nx,1,1,...] - ALLOCATE (ztile, int, header->Naxes); - if (!gfits_scan (header, "ZTILE1", "%d", 1, &ztile[0])) { - ztile[0] = header->Naxis[0]; - for (i = 1; i < header->Naxes; i++) { - ztile[i] = 1; - } - } else { - gfits_delete (header, "ZTILE1", 1); - for (i = 1; i < header->Naxes; i++) { - snprintf (key, 10, "ZTILE%d", i + 1); - if (!gfits_scan (header, key, "%d", 1, &ztile[i])) ESCAPE; - gfits_delete (header, key, 1); - } - } + snprintf (keyword, 10, "ZTILE%d", i + 1); + 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: + // ZEXTEND + // ZBLOCKED + // ZHECKSUM + // ZDATASUM + // PCOUNT -> ZPCOUNT + // GCOUNT -> ZGCOUNT + + float zscale = 1.0; + // gfits_modify (header, "ZSCALE", "%f", 1, zscale); + + float zzero = 0.0; + // gfits_modify (header, "ZZERO", "%f", 1, zzero); + + int zblank = 32767; + // gfits_modify (header, "ZBLANK", "%d", 1, zblank); + + int oblank = 32767; + // gfits_modify (header, "BLANK", "%d", 1, oblank); + +# if (0) + // need to set ZNAMEnn and ZVALnn if they are defined and used // search for algorithm-specific keywords. these are used to control compression options // note the difference between the keywords (1 indexed) and the variables (0 indexed) - NOPTIONS = 10; - ALLOCATE (optname, char *, NOPTIONS); - ALLOCATE (optvalue, char *, NOPTIONS); - for (Noptions = 0; TRUE; Noptions++) { - snprintf (key, 10, "ZNAME%d", Noptions + 1); - if (!gfits_scan (header, key, "%s", 1, word)) break; - gfits_delete (header, key, 1); - optname[Noptions] = strcreate (word); - - snprintf (key, 10, "ZVAL%d", Noptions + 1); - if (!gfits_scan (header, key, "%s", 1, word)) ESCAPE; - gfits_delete (header, key, 1); - optvalue[Noptions] = strcreate (word); - - if (Noptions == NOPTIONS - 1) { - NOPTIONS += 10; - REALLOCATE (optname, char *, NOPTIONS); - REALLOCATE (optvalue, char *, NOPTIONS); - } - } - - // XXX get ZMASKCMP - int zsimple, have_zsimple; - int have_ztension; - - // 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 - if (have_zsimple && have_ztension) ESCAPE; - - // if neither are present, we have an image that is not really following the standard: - // assume it is a PHU - if (!have_zsimple && !have_ztension) { - header->simple = TRUE; - gfits_extended_to_primary (header, header->simple, "Image data"); - - MOD_KEYWORD_ALT ("ZEXTEND", "EXTEND", "%t", &header->extend, header->extend); - MOD_KEYWORD_ALT ("ZBLOCKED", "BLOCKED", "%t", &header->extend, header->extend); - } - - // have_zsimple : image comes from a PHU - if (have_zsimple) { - header->simple = zsimple; - gfits_delete (header, "ZSIMPLE", 1); - gfits_extended_to_primary (header, header->simple, "Image data"); - - MOD_KEYWORD_ALT ("ZEXTEND", "EXTEND", "%t", &header->extend, header->extend); - MOD_KEYWORD_ALT ("ZBLOCKED", "BLOCKED", "%t", &header->extend, header->extend); - } - - // have_ztension : image comes from an extension - if (have_ztension) { - gfits_delete (header, "ZTENSION", 1); - gfits_modify_extended (header, exttype, "Image extension"); - - // we may have an uncompressed PCOUNT / GCOUNT value, otherwise set to 0,1 - if (gfits_scan (header, "ZPCOUNT", "%d", 1, &header->pcount)) { - gfits_delete (header, "ZPCOUNT", 1); - gfits_modify (header, "PCOUNT", "%d", 1, header->pcount); + for (i = 0; i < Noptions; i++) { + snprintf (key, 10, "ZNAME%d", i + 1); + if (!gfits_modify (header, key, "%s", 1, optname[i])) break; + + snprintf (key, 10, "ZVAL%d", i + 1); + if (!gfits_modify (header, key, "%s", 1, optvalue[i])) break; + } +# endif + + VarLengthColumn zdef; + 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 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); + + // define compression-specific keywords, update header as needed. + if (!strcasecmp(zcmptype, "RICE_1") || !strcasecmp(zcmptype, "RICE_ONE")) { + 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 +# if VERBOSE_DUMP + int tile_pixsize = abs(header[0].bitpix) / 8; + fprintf (stderr, "raw_pixsize: %d, cmp_pixsize: %d, tile_pixsize: %d, raw_bitpix: %d\n", raw_pixsize, cmp_pixsize, tile_pixsize, raw_bitpix); +# endif + + // allocate the buffer for compression work + int Nzdata_alloc = raw_pixsize*max_tile_size + 100; + ALLOCATE (raw, char, raw_pixsize*max_tile_size); + ALLOCATE (zdata, char, Nzdata_alloc); + + // init the otile[] counters + if (!gfits_imtile_start (matrix, otile)) ESCAPE; + + // compress the data : copy into a tile, compress the tile, then add to the output table + // each tile -> 1 row of the output table + for (i = 0; i < Ntile; i++) { + + // 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 + // for float -> int scaling by zscale, zzero may be applied + if (!strcasecmp (zcmptype, "GZIP_2") || !strcasecmp (zcmptype, "NONE_2")) { + if (!gfits_collect_gzp2 (matrix, raw, Nraw, raw_bitpix, otile, oblank, ztile, zblank, zscale, zzero)) ESCAPE; } else { - header->pcount = 0; - gfits_modify (header, "PCOUNT", "%d", 1, header->pcount); - } - if (gfits_scan (header, "ZGCOUNT", "%d", 1, &header->gcount)) { - gfits_delete (header, "ZGCOUNT", 1); - gfits_modify (header, "GCOUNT", "%d", 1, header->gcount); - } else { - header->pcount = 1; - gfits_modify (header, "GCOUNT", "%d", 1, header->gcount); - } - } else { - header->pcount = 0; - header->gcount = 1; - gfits_modify (header, "PCOUNT", "%d", 1, header->pcount); - gfits_modify (header, "GCOUNT", "%d", 1, header->gcount); - } - - MOD_KEYWORD ("ZHECKSUM", "CHECKSUM", "%s", checksum, checksum); - 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); - - // 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, word)) ESCAPE; - if (!strcmp (word, "COMPRESSED_DATA")) break; - } - zcol = i; - - if (!gfits_varlength_column_define (ftable, &zdef, zcol)) ESCAPE; - gfits_delete (header, "TFIELDS", 1); - snprintf (key, 10, "TTYPE"OFF_T_FMT, zcol); - gfits_delete (header, key, 1); - snprintf (key, 10, "TFORM"OFF_T_FMT, zcol); - gfits_delete (header, key, 1); - - // create the output image - gfits_create_matrix (header, matrix); - - // 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]; - } - - // this takes place in three major steps: - // 1) read the table data : zdef.format tells the size of the varlength heap element - // 2) uncompress the data : pixel size depends on compression method - // 3) distribute data to the tiles : output pixel size depends on header.bitpix - - // heapdata -> zdata -> odata -> idata - - // heapdata.pixsize : zdef.format - // zdata.pixsize : depends on compression - // odata.pixsize : depends on compression - // idata.pixsize : header.bitpix - - // size of an element in the vartable heap section - zdata_pixsize = gfits_vartable_heap_pixsize (zdef.format); - - // size of a pixel in the final image (not needed) - // idata_pixsize = abs(header[0].bitpix) / 8; - - // size of a pixel in the output from the decompression routine - odata_pixsize = gfits_uncompressed_data_pixsize (cmptype, header[0].bitpix, optname, optvalue, Noptions); - ALLOCATE (out, char, odata_pixsize*max_tile_size); - - off_t row; - - // uncompress the data - Nzrows = ftable->header->Naxis[1]; - for (row = 0; row < Nzrows; row++) { - - // expected output size for this tile - Nout = gfits_tile_size (matrix, otile, ztile); - - zdata = gfits_varlength_column_pointer (ftable, &zdef, row, &Nzdata); - if (!zdata) return (FALSE); - - // 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)) return (FALSE); - } - - // 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)) return (FALSE); - - // copy the uncompressed pixels into their correct locations - if (!gfits_distribute_data (matrix, odata_pixsize, out, Nout, otile, oblank, ztile, zblank, zscale, zzero)) return (FALSE); - - // 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; - } - } - } + if (!gfits_collect_data (matrix, raw, Nraw, raw_bitpix, otile, oblank, ztile, zblank, zscale, zzero)) ESCAPE; + } + + if (VERBOSE_DUMP && (i == 0)) { + int k; + fprintf (stderr, "cmp mat: "); + for (k = 0; k < 64; k++) { fprintf (stderr, "%02hhx", matrix->buffer[k]); if (k % 4 == 3) fprintf (stderr, " "); } + fprintf (stderr, "\n"); + fprintf (stderr, "cmp raw: "); + for (k = 0; k < 64; k++) { fprintf (stderr, "%02hhx", raw[k]); if (k % 4 == 3) fprintf (stderr, " "); } + fprintf (stderr, "\n"); + } + + // gzip compresses bytes which are in BIG-ENDIAN order + // Nraw is number of pixels + if (!strcasecmp(zcmptype, "GZIP_1")) { + if (!gfits_byteswap_zdata (raw, Nraw * raw_pixsize, raw_pixsize)) ESCAPE; + } + int Nzdata = Nzdata_alloc; // available space, replaced with actual output size on compression + + if (VERBOSE_DUMP && (i == 0)) { + int k; + fprintf (stderr, "cmp swp: "); + for (k = 0; k < 64; k++) { fprintf (stderr, "%02hhx", raw[k]); if (k % 4 == 3) fprintf (stderr, " "); } + fprintf (stderr, "\n"); + fprintf (stderr, "Nzdata: %d -> ", Nzdata); + } + + // optname, optvalue = NULL, Noptions = 0 : defaults for RICE_1, HCOMPRESS_1 + + if (!gfits_compress_data (zdata, &Nzdata, zcmptype, NULL, NULL, 0, raw, Nraw, raw_pixsize, ztile[1], ztile[0])) ESCAPE; + if (VERBOSE_DUMP && (i == 0)) { + int k; + fprintf (stderr, "%d\n", Nzdata); + fprintf (stderr, "cmp SWP: "); + for (k = 0; k < 64; k++) { fprintf (stderr, "%02hhx", zdata[k]); if (k % 4 == 3) fprintf (stderr, " "); } + fprintf (stderr, "\n"); + } + + // 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, "NONE_2") && // strcasecmp(zcmptype, "NONE") && + strcasecmp(zcmptype, "GZIP_1") && + strcasecmp(zcmptype, "GZIP_2") && + strcasecmp(zcmptype, "RICE_1") && + strcasecmp(zcmptype, "RICE_ONE")) { + // Nzdata is number of bytes + if (!gfits_byteswap_zdata (zdata, Nzdata, cmp_pixsize)) ESCAPE; + } + + if (VERBOSE_DUMP && (i == 0)) { + int k; + fprintf (stderr, "cmp dat: "); + for (k = 0; k < 64; k++) { fprintf (stderr, "%02hhx", zdata[k]); if (k % 4 == 3) fprintf (stderr, " "); } + fprintf (stderr, "\n"); + } + if (!gfits_varlength_column_add_data (ftable, zdata, Nzdata, i, &zdef)) ESCAPE; + + // update the otile[] counters, carrying to the next dimension if needed + gfits_imtile_next (matrix, ztile, ntile, otile); + } + + if (!gfits_varlength_column_finish (ftable, &zdef)) ESCAPE; + if (VERBOSE_DUMP) { + int k; + fprintf (stderr, "cmp tbl: "); + for (k = 0; k < 64; k++) { fprintf (stderr, "%02hhx", ftable->buffer[k]); if (k % 4 == 3) fprintf (stderr, " "); } + fprintf (stderr, "\n"); + } + + free (raw); + free (zdata); + + free (otile); + free (ntile); + + if (!Ztile) free (ztile); return (TRUE); + +escape: + FREE (raw); + FREE (zdata); + + FREE (otile); + FREE (ntile); + + if (!Ztile) FREE (ztile); + return FALSE; } -// gfits_varlength returns a pointer to a chunk of Nzdata of data elements starting at zdata. - -// 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. - -// XXX inconsistency between data element size between compressed data, uncompressed data, -// and output pixel data..... - -// true sizes of this tile (in pixels) -off_t gfits_tile_size (Matrix *matrix, int *otile, int *ztile) { - - off_t i, Npixels, Ndimen; - - Npixels = 1; - for (i = 0; i < matrix->Naxes; i++) { - Ndimen = MIN ((matrix->Naxis[i] - otile[i]*ztile[i]), ztile[i]); - Npixels *= Ndimen; - } - - return (Npixels); +int gfits_copy_keywords_compress (Header *srchead, Header *tgthead) { + + int i; + + for (i = 0; i < srchead->datasize; i += FT_LINE_LENGTH) { + } + return TRUE; } -// 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) { +// 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_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; @@ -316,9 +309,7 @@ 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) + 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; @@ -332,8 +323,8 @@ } - // double check reported size (needs the Isize) - assert (Ndata == Ztile[0]*Nline); - - // set the starting point of the tile: + // 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] + ...)); @@ -347,52 +338,82 @@ offset = 0; + int directCopy = (zzero == 0.0) && (zscale == 1.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++) { \ - if (*Iptr == zblank) { \ - *Optr = oblank; \ - } else { \ - *Optr = *Iptr * zscale + zzero; \ +# define SCALE_AND_DIST_INT_PRINT(TYPE, SIZE, INTYPE) { \ + TYPE *TILEptr = (TYPE *) &matrix->buffer[SIZE*(offset + start)]; \ + for (j = 0; j < Ztile[0]; j++, TILEptr++, RAWptr++) { \ + if (*TILEptr == oblank) { \ + *RAWptr = zblank; \ + } else { \ + /* *RAWptr = directCopy ? *TILEptr : (*TILEptr - zzero) / zscale; */ \ + *RAWptr = *TILEptr; \ + fprintf (stderr, "col: %s to %s : %08x : %08x\n", #TYPE, #INTYPE, *TILEptr, *RAWptr); \ } } } - // 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(); \ + // 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_INT(OTYPE, OSIZE) { \ + OTYPE *TILEptr = (OTYPE *) &matrix->buffer[OSIZE*(offset + start)]; \ + for (j = 0; j < Ztile[0]; j++, TILEptr++, RAWptr++) { \ + if (*TILEptr == oblank) { \ + *RAWptr = zblank; \ + } else { \ + if (directCopy) { \ + *RAWptr = *TILEptr; \ + } else { \ + *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 +# define SCALE_AND_DIST_FLOAT(OTYPE, OSIZE) { \ + OTYPE *TILEptr = (OTYPE *) &matrix->buffer[OSIZE*(offset + start)]; \ + for (j = 0; j < Ztile[0]; j++, TILEptr++, RAWptr++) { \ + if (!isfinite(*TILEptr)) { \ + *RAWptr = zblank; \ + } else { \ + if (directCopy) { \ + *RAWptr = *TILEptr; \ + } else { \ + *RAWptr = (*TILEptr - zzero) / zscale; \ + } } } } + + // this macro sets up the outer switch and calls above the macro with all RAW buffer bitpix options + // +# define SETUP_RAWSIZE_PRINT(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_PRINT (int, 4, TYPE); break; \ + case -32: SCALE_AND_DIST_FLOAT (float, 4); break; \ + case -64: SCALE_AND_DIST_FLOAT (double, 8); break; \ + default: abort(); \ + } } + + // this macro sets up the outer switch and calls above the macro with all RAW buffer bitpix options + // +# define SETUP_RAWSIZE(ITYPE, ISIZE) { \ + ITYPE *RAWptr = (ITYPE *) &raw[i*ISIZE*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; - 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(); } @@ -421,123 +442,107 @@ free (counter); free (Ztile); + return (TRUE); } -int gfits_byteswap_zdata (char *zdata, int Nzdata, int pixsize) { - -# define DOSWAP(A,B) { char tmp = A; A = B; B = tmp; } - -# ifdef BYTE_SWAP - - int i; - - switch (pixsize) { - case 1: - break; - - case 2: - for (i = 0; i < 2*Nzdata; i+=2) { - DOSWAP (zdata[i+0], zdata[i+1]); +// 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*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"); } - break; - - case 4: - for (i = 0; i < 4*Nzdata; i+=4) { - DOSWAP (zdata[i+0], zdata[i+3]); - DOSWAP (zdata[i+1], zdata[i+2]); + + // 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; + } } - break; - - case 8: - for (i = 0; i < 8*Nzdata; i+=8) { - DOSWAP (zdata[i+0], zdata[i+7]); - DOSWAP (zdata[i+1], zdata[i+6]); - DOSWAP (zdata[i+2], zdata[i+5]); - DOSWAP (zdata[i+3], zdata[i+4]); - } - break; - } -# endif + 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); } - -int gfits_extension_is_compressed (Header *header) { - - int has_extension, has_extname, has_zimage, zimage; - char extname[80], extension[80]; - - has_extension = gfits_scan (header, "XTENSION", "%s", 1, extension); - has_extname = gfits_scan (header, "EXTNAME", "%s", 1, extname); - has_zimage = gfits_scan_alt (header, "ZIMAGE", "%t", 1, &zimage); - - if (has_extension && !strcmp (extension, "IMAGE")) return (FALSE); - 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 zimage = FALSE; - // int ztension = FALSE; - - int has_zimage = gfits_scan_alt (header, "ZIMAGE", "%t", 1, &zimage); - // int has_ztension = gfits_scan_alt (header, "ZTENSION", "%t", 1, &ztension); - - if (has_zimage && zimage) return (TRUE); - - return (FALSE); -} - -int gfits_uncompressed_data_pixsize (char *cmptype, int out_bitpix, char **optname, char **optvalue, int Noptions) { - - int i, Nbyte; - - if (!strcasecmp(cmptype, "GZIP_1")) { - 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")) { - - // 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); - } - } - - return (4); - } - 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_vartable_heap_pixsize (char format) { - - if (format == 'B') { - return (1); - } - if (format == 'I') { - return (2); - } - if (format == 'J') { - return (4); - } - fprintf (stderr, "invalid size for compressed data: %c\n", format); - abort (); -} Index: trunk/Ohana/src/libfits/matrix/F_compress_data.c =================================================================== --- trunk/Ohana/src/libfits/matrix/F_compress_data.c (revision 38441) +++ trunk/Ohana/src/libfits/matrix/F_compress_data.c (revision 38441) @@ -0,0 +1,159 @@ +# include +# include +# include + +/* 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); +int fits_rdecomp_byte (unsigned char *c, int clen, unsigned char array[], int nx, int nblock); + +/* functions defined in fits_hcompress.c */ +# define LONGLONG long long +int fits_hcompress(int *a, int ny, int nx, int scale, char *output, long *nbytes, int *status); +int fits_hcompress64(LONGLONG *a, int ny, int nx, int scale, char *output, long *nbytes, int *status); + +/* functions defined in fits_hdeccompress.c */ +int fits_hdecompress(unsigned char *input, int smooth, int *a, int *ny, int *nx, int *scale, int *status); +int fits_hdecompress64(unsigned char *input, int smooth, LONGLONG *a, int *ny, int *nx, int *scale, int *status); + +/* functions defined in pliocomp.c */ +int pl_p2li (int *pxsrc, int xs, short *lldst, int npix); +int pl_l2pi (short *ll_src, int xs, int *px_dst, int npix); + +/* functions defined in gzip.c */ +int gfits_gz_stripheader (unsigned char *data, int *ndata); +int gfits_uncompress (Bytef *dest, uLongf *destLen, const Bytef *source, uLong sourceLen); +int gfits_compress (Bytef *dest, uLongf *destLen, const Bytef *source, uLong sourceLen); + +# 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, + int Nx, int Ny) { + + int Nout; + int status; + + // do not actually compress : this is used for testing + if (!strcasecmp(cmptype, "NONE") || !strcasecmp(cmptype, "NONE_1") || !strcasecmp(cmptype, "NONE_2")) { + unsigned long Nbytes = Nrawpix * rawpix_size; + + memcpy (zdata, rawdata, Nbytes); + *Nzdata = Nbytes; + return (TRUE); + } + + // 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; + + // the data must be byteswapped before compression begins + // if (!gfits_byteswap_zdata (rawdata, Nrawpix, rawpix_size)) return (FALSE); + + // Nzdata is the size of the allocated buffer before the function is called, returns the actual size + // Nbytes is the number of bytes in the raw data buffer + uLongf destLen = *Nzdata; + status = gfits_compress ((Bytef *) zdata, &destLen, (Bytef *) rawdata, Nbytes); + *Nzdata = destLen; + if (status != Z_OK) ESCAPE(FALSE); + return (TRUE); + } + + if (!strcasecmp(cmptype, "RICE_1") || !strcasecmp(cmptype, "RICE_ONE")) { + int i, blocksize; + // look for the BLOCKSIZE + blocksize = 32; + for (i = 0; i < Nopt; i++) { + if (!strcmp(optname[i], "BLOCKSIZE")) { + blocksize = atoi (optvalue[i]); + if ((blocksize != 16) && (blocksize != 32)) { + fprintf (stderr, "RICE blocksize is not valid: %d (%s = %s)\n", blocksize, optname[i], optvalue[i]); + return (FALSE); + } + } + } + + // Ninsum += Nzdata; + // Noutsum += *Nout; + // fprintf (stderr, "%d comp bytes; %d uncomp 'pixels', totals: %d %d\n", Nzdata, Npix, Ninsum, Noutsum); + + int k; + + switch (rawpix_size) { + case 4: + // rice compression from the CFITSIO source tree + Nout = fits_rcomp ((int *) rawdata, Nrawpix, (unsigned char *) zdata, *Nzdata, blocksize); + break; + + case 2: + // Nout is the number of bytes in the compressed buffer + Nout = fits_rcomp_short ((short *) rawdata, Nrawpix, (unsigned char *) zdata, *Nzdata, blocksize); + + if (0) { + fprintf (stderr, "Nout: %d, Nrawpix: %d\n", Nout, Nrawpix); + fprintf (stderr, "cmp inp: "); + for (k = 0; k < Nout; k++) { fprintf (stderr, "%02hhx", zdata[k]); if (k % 2) fprintf (stderr, " "); } + fprintf (stderr, "\n"); + } + break; + + case 1: + Nout = fits_rcomp_byte ((char *) rawdata, Nrawpix, (unsigned char *) zdata, *Nzdata, blocksize); + break; + + default: + 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")) { + // 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")) { + long Nbytes = Nrawpix * rawpix_size; + int scale = 0; + status = 0; + // call hdecompress without smoothing + + 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"); + return (FALSE); + } + // fprintf (stderr, "decompression yields image %d x %d (scale: %d)\n", Nx, Ny, scale); + + *Nzdata = Nbytes / rawpix_size; + return (TRUE); + } + + fprintf (stderr, "unknown compression %s\n", cmptype); + return (FALSE); +} + Index: trunk/Ohana/src/libfits/matrix/F_compress_utils.c =================================================================== --- trunk/Ohana/src/libfits/matrix/F_compress_utils.c (revision 38441) +++ trunk/Ohana/src/libfits/matrix/F_compress_utils.c (revision 38441) @@ -0,0 +1,344 @@ +# include +# include +# include + +// gfits_varlength returns a pointer to a chunk of Nzdata of data elements starting at zdata. + +// 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. + +// XXX inconsistency between data element size between compressed data, uncompressed data, +// and output pixel data..... + +// need a structure to describe tiles, current tile? + +typedef struct { + int *ztile; // max size of a tile (edge tiles may be smaller) + int *otile; // counter for current tile (eg, for (2,3,0) otile[0] = 2, otile[1] = 3, otile[2] = 0) + int *ntile; // number of tiles in dimension [i] +} ImageTile; + +// advance the otile counter. +int gfits_imtile_next (Matrix *matrix, int *ztile, int *ntile, int *otile) { + + int i; + + // update the tile counters, carrying to the next dimension if needed + for (i = 0; i < matrix->Naxes; i++) { + otile[i] ++; + if (otile[i] == ntile[i]) { + if (i == matrix->Naxes - 1) return FALSE; + otile[i] = 0; + } else { + return TRUE; + } + } + return TRUE; +} + +// how many tiles needed for this image (given ztile[]) +int gfits_imtile_start (Matrix *matrix, int *otile) { + + int i; + + // update the tile counters, carrying to the next dimension if needed + for (i = 0; i < matrix->Naxes; i++) { + otile[i] = 0; + } + return TRUE; +} + +// how many tiles needed for this image (given ztile[]) +int gfits_imtile_count (Matrix *matrix, int *ztile, int *ntile) { + + int i; + + int Ntiles = 1; + for (i = 0; i < matrix->Naxes; i++) { + ntile[i] = (matrix->Naxis[i] % ztile[i]) ? (matrix->Naxis[i] / ztile[i] + 1) : (matrix->Naxis[i] / ztile[i]); + Ntiles *= ntile[i]; + } + + return (Ntiles); +} + +// how many tiles needed for this image (given ztile[]) +int gfits_imtile_maxsize (Matrix *matrix, int *ztile) { + + int i; + + int max_tile_size = 1; + for (i = 0; i < matrix->Naxes; i++) { + max_tile_size *= ztile[i]; + } + + return (max_tile_size); +} + +// true sizes of this tile (in pixels) +off_t gfits_tile_size (Matrix *matrix, int *otile, int *ztile) { + + off_t i, Npixels, Ndimen; + + Npixels = 1; + for (i = 0; i < matrix->Naxes; i++) { + Ndimen = MIN ((matrix->Naxis[i] - otile[i]*ztile[i]), ztile[i]); + Npixels *= Ndimen; + } + + return (Npixels); +} + +int gfits_byteswap_zdata (char *zdata, int Nzdata, int pixsize) { + +# define DOSWAP(A,B) { char tmp = A; A = B; B = tmp; } + +# ifdef BYTE_SWAP + + int i; + + // fprintf (stderr, "swapping %d bytes in pix of size %d bytes...\n", Nzdata, pixsize); + switch (pixsize) { + case 1: + break; + + case 2: + for (i = 0; i < Nzdata; i+=2) { + DOSWAP (zdata[i+0], zdata[i+1]); + } + break; + + case 4: + for (i = 0; i < Nzdata; i+=4) { + DOSWAP (zdata[i+0], zdata[i+3]); + DOSWAP (zdata[i+1], zdata[i+2]); + } + break; + + case 8: + for (i = 0; i < Nzdata; i+=8) { + DOSWAP (zdata[i+0], zdata[i+7]); + DOSWAP (zdata[i+1], zdata[i+6]); + DOSWAP (zdata[i+2], zdata[i+5]); + DOSWAP (zdata[i+3], zdata[i+4]); + } + break; + } +# endif + return (TRUE); +} + +int gfits_extension_is_compressed_image (Header *header) { + + int has_extension, has_extname, has_zimage, zimage; + char extname[80], extension[80]; + + has_extension = gfits_scan (header, "XTENSION", "%s", 1, extension); + has_extname = gfits_scan (header, "EXTNAME", "%s", 1, extname); + has_zimage = gfits_scan_alt (header, "ZIMAGE", "%t", 1, &zimage); + + if (has_extension && !strcmp (extension, "IMAGE")) return (FALSE); + if (has_zimage && zimage) return (TRUE); + if (has_extname) { + if (!strcmp (extname, "COMPRESSED_IMAGE")) return (TRUE); + } + + return (FALSE); +} + +int gfits_extension_is_compressed_table (Header *header) { + + int has_ztable, ztable; + has_ztable = gfits_scan_alt (header, "ZTABLE", "%t", 1, &ztable); + if (has_ztable && ztable) return (TRUE); + return (FALSE); +} + +int gfits_compressed_is_primary (Header *header) { + + int zimage = FALSE; + // int ztension = FALSE; + + int has_zimage = gfits_scan_alt (header, "ZIMAGE", "%t", 1, &zimage); + // int has_ztension = gfits_scan_alt (header, "ZTENSION", "%t", 1, &ztension); + + if (has_zimage && zimage) return (TRUE); + + return (FALSE); +} + +int gfits_cmptype_valid (char *cmptype) { + + if (!strcasecmp(cmptype, "NONE")) return TRUE; + if (!strcasecmp(cmptype, "NONE_1")) return TRUE; // do not compress, but shuffle a la GZIP_1 + if (!strcasecmp(cmptype, "NONE_2")) return TRUE; // do not compress, but shuffle a la GZIP_2 + if (!strcasecmp(cmptype, "GZIP_1")) return TRUE; + if (!strcasecmp(cmptype, "GZIP_2")) return TRUE; + if (!strcasecmp(cmptype, "PLIO_1")) return TRUE; + if (!strcasecmp(cmptype, "RICE_1")) return TRUE; + if (!strcasecmp(cmptype, "RICE_ONE")) return TRUE; + if (!strcasecmp(cmptype, "HCOMPRESS_1")) return TRUE; + return FALSE; +} + +int gfits_uncompressed_data_pixsize (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_2") || + !strcasecmp(cmptype, "NONE_1") || + !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); + } + + // 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_compressed_data_pixsize (char *cmptype, int out_bitpix, char **optname, char **optvalue, int Noptions) { + + int i, Nbyte; + + // NONE should be swapped to mimic swapping types + if (!strcasecmp(cmptype, "NONE") || !strcasecmp(cmptype, "NONE_1")) { + 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); + } + + // GZIP_1, GZIP_2 should not be swapped after compression + if (!strcasecmp(cmptype, "GZIP_1") || + !strcasecmp(cmptype, "GZIP_2") || + !strcasecmp(cmptype, "NONE_2")) { + 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_1") || + !strcasecmp(cmptype, "NONE_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) { + + if (format == 'B') { + return (1); + } + if (format == 'I') { + return (2); + } + if (format == 'J') { + return (4); + } + fprintf (stderr, "invalid size for compressed data: %c\n", format); + abort (); +} Index: trunk/Ohana/src/libfits/matrix/F_matrix.c =================================================================== --- trunk/Ohana/src/libfits/matrix/F_matrix.c (revision 38255) +++ trunk/Ohana/src/libfits/matrix/F_matrix.c (revision 38441) @@ -2,29 +2,5 @@ # include -off_t gfits_data_size (Header *header) { - - int i; - off_t Nrec, size; - - if (header[0].Naxes == 0) return (0); - - size = abs(header[0].bitpix / 8); - - for (i = 0; i < header[0].Naxes; i++) - size *= header[0].Naxis[i]; - - // XXX do I multiply this times gcount? - size += header[0].pcount; - - /* round up to next complete block */ - if (size % FT_RECORD_SIZE) { - Nrec = 1 + (int) (size / FT_RECORD_SIZE); - size = FT_RECORD_SIZE * Nrec; - } - - return (size); -} - -off_t gfits_data_min_size (Header *header) { +off_t gfits_heap_start (Header *header) { int i; @@ -38,6 +14,13 @@ size *= header[0].Naxis[i]; - // XXX do I multiply this times gcount? + return (size); +} + +off_t gfits_data_min_size (Header *header) { + + off_t size = gfits_heap_start (header); size += header[0].pcount; + + // XXX what do I do with gcount? return (size); @@ -46,11 +29,9 @@ off_t gfits_data_pad_size (off_t rawsize) { - off_t Nrec, size; - - size = rawsize; + off_t size = rawsize; /* round up to next complete block */ if (rawsize % FT_RECORD_SIZE) { - Nrec = 1 + (int) (rawsize / FT_RECORD_SIZE); + off_t Nrec = 1 + (int) (rawsize / FT_RECORD_SIZE); size = FT_RECORD_SIZE * Nrec; } @@ -58,2 +39,10 @@ return (size); } + +off_t gfits_data_size (Header *header) { + off_t rawsize = gfits_data_min_size(header); + off_t size = gfits_data_pad_size (rawsize); + + return (size); +} + Index: trunk/Ohana/src/libfits/matrix/F_uncompress_M.c =================================================================== --- trunk/Ohana/src/libfits/matrix/F_uncompress_M.c (revision 38441) +++ trunk/Ohana/src/libfits/matrix/F_uncompress_M.c (revision 38441) @@ -0,0 +1,603 @@ +# include +# include +# 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 { \ + fprintf (stderr, "error in %s @ line %d\n", __func__, __LINE__); \ + 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 (FALSE); } + +# define MOD_KEYWORD(ZNAME,NAME,TYPE,IN,OUT) { \ + if (gfits_scan (header, ZNAME, TYPE, 1, IN)) { \ + gfits_modify (header, NAME, TYPE, 1, OUT); \ + } \ + gfits_delete (header, ZNAME, 1); } + +# define MOD_KEYWORD_ALT(ZNAME,NAME,TYPE,IN,OUT) { \ + if (gfits_scan_alt (header, ZNAME, TYPE, 1, IN)) { \ + gfits_modify_alt (header, NAME, TYPE, 1, OUT); \ + } \ + gfits_delete (header, ZNAME, 1); } + +# define MOD_KEYWORD_REQUIRED(ZNAME,NAME,TYPE,IN,OUT) { \ + if (!gfits_scan (header, ZNAME, TYPE, 1, IN)) ESCAPE; \ + gfits_delete (header, ZNAME, 1); \ + gfits_modify (header, NAME, TYPE, 1, OUT); } + +int gfits_uncompress_image (Header *header, Matrix *matrix, FTable *ftable) { + + off_t Nzdata, Nzrows, zcol; + int i, j, status, zimage, Nout, max_tile_size; + char cmptype[80]; + char zaxis[10], naxis[10], key[10], word[81], exttype[81], checksum[81], datasum[81]; + int Noptions, NOPTIONS, zblank, oblank; + VarLengthColumn zdef; + float zscale, zzero; + + int *ztile = NULL; + int *otile = NULL; + int *ntile = NULL; + char **optname = NULL; + char **optvalue = NULL; + char *out = NULL; + char *zdata = NULL; + + Noptions = 0; + + // is ZIMAGE present? + // NOTE target of %t must be int length + status = gfits_scan_alt (ftable->header, "ZIMAGE", "%t", 1, &zimage); + if (!status || !zimage) ESCAPE; + + // copy original header to output header + gfits_copy_header (ftable->header, header); + + // delete ZIMAGE from output header + gfits_delete (header, "ZIMAGE", 1); + + // extract compression-specific keywords, update header as needed. + if (!gfits_scan (header, "ZCMPTYPE", "%s", 1, cmptype)) ESCAPE; + gfits_delete (header, "ZCMPTYPE", 1); + + MOD_KEYWORD_REQUIRED ("ZBITPIX", "BITPIX", "%d", &header->bitpix, header->bitpix); + MOD_KEYWORD_REQUIRED ("ZNAXIS", "NAXIS", "%d", &header->Naxes, header->Naxes); + + for (i = 0; i < header->Naxes; i++) { + snprintf (zaxis, 10, "ZNAXIS%d", i + 1); + snprintf (naxis, 10, "NAXIS%d", i + 1); + MOD_KEYWORD_REQUIRED (zaxis, naxis, OFF_T_FMT, &header->Naxis[i], header->Naxis[i]); + } + + // set up the tile sizes : ztile is the default size of the tile in the Nth dimension + // the actual tile size may be smaller at the edge of a dimension. if the ZTILEn + // entries are not found, default to [Nx,1,1,...] + ALLOCATE (ztile, int, header->Naxes); + if (!gfits_scan (header, "ZTILE1", "%d", 1, &ztile[0])) { + ztile[0] = header->Naxis[0]; + for (i = 1; i < header->Naxes; i++) { + ztile[i] = 1; + } + } else { + gfits_delete (header, "ZTILE1", 1); + for (i = 1; i < header->Naxes; i++) { + snprintf (key, 10, "ZTILE%d", i + 1); + if (!gfits_scan (header, key, "%d", 1, &ztile[i])) ESCAPE; + gfits_delete (header, key, 1); + } + } + + // search for algorithm-specific keywords. these are used to control compression options + // note the difference between the keywords (1 indexed) and the variables (0 indexed) + NOPTIONS = 10; + ALLOCATE (optname, char *, NOPTIONS); + ALLOCATE (optvalue, char *, NOPTIONS); + for (Noptions = 0; TRUE; Noptions++) { + snprintf (key, 10, "ZNAME%d", Noptions + 1); + if (!gfits_scan (header, key, "%s", 1, word)) break; + gfits_delete (header, key, 1); + optname[Noptions] = strcreate (word); + + snprintf (key, 10, "ZVAL%d", Noptions + 1); + if (!gfits_scan (header, key, "%s", 1, word)) ESCAPE; + gfits_delete (header, key, 1); + optvalue[Noptions] = strcreate (word); + + if (Noptions == NOPTIONS - 1) { + NOPTIONS += 10; + REALLOCATE (optname, char *, NOPTIONS); + REALLOCATE (optvalue, char *, NOPTIONS); + } + } + + // XXX get ZMASKCMP + int zsimple, have_zsimple; + int have_ztension; + + // 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 + if (have_zsimple && have_ztension) ESCAPE; + + // if neither are present, we have an image that is not really following the standard: + // assume it is a PHU + if (!have_zsimple && !have_ztension) { + header->simple = TRUE; + gfits_extended_to_primary (header, header->simple, "Image data"); + + MOD_KEYWORD_ALT ("ZEXTEND", "EXTEND", "%t", &header->extend, header->extend); + MOD_KEYWORD_ALT ("ZBLOCKED", "BLOCKED", "%t", &header->extend, header->extend); + } + + // have_zsimple : image comes from a PHU + if (have_zsimple) { + header->simple = zsimple; + gfits_delete (header, "ZSIMPLE", 1); + gfits_extended_to_primary (header, header->simple, "Image data"); + + MOD_KEYWORD_ALT ("ZEXTEND", "EXTEND", "%t", &header->extend, header->extend); + MOD_KEYWORD_ALT ("ZBLOCKED", "BLOCKED", "%t", &header->extend, header->extend); + } + + // have_ztension : image comes from an extension + if (have_ztension) { + gfits_delete (header, "ZTENSION", 1); + gfits_modify_extended (header, exttype, "Image extension"); + + // we may have an uncompressed PCOUNT / GCOUNT value, otherwise set to 0,1 + if (gfits_scan (header, "ZPCOUNT", "%d", 1, &header->pcount)) { + gfits_delete (header, "ZPCOUNT", 1); + gfits_modify (header, "PCOUNT", "%d", 1, header->pcount); + } else { + header->pcount = 0; + gfits_modify (header, "PCOUNT", "%d", 1, header->pcount); + } + if (gfits_scan (header, "ZGCOUNT", "%d", 1, &header->gcount)) { + gfits_delete (header, "ZGCOUNT", 1); + gfits_modify (header, "GCOUNT", "%d", 1, header->gcount); + } else { + header->pcount = 1; + gfits_modify (header, "GCOUNT", "%d", 1, header->gcount); + } + } else { + header->pcount = 0; + header->gcount = 1; + gfits_modify (header, "PCOUNT", "%d", 1, header->pcount); + gfits_modify (header, "GCOUNT", "%d", 1, header->gcount); + } + + MOD_KEYWORD ("ZHECKSUM", "CHECKSUM", "%s", checksum, checksum); + 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); + + // find the COMPRESSED_DATA column (format should be 1PB, 1PI, 1PJ) + // is it required that this be the only column? + for (i = 1; TRUE; i++) { + snprintf (key, 10, "TTYPE%d", i); + 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; + gfits_delete (header, "TFIELDS", 1); + snprintf (key, 10, "TTYPE"OFF_T_FMT, zcol); + gfits_delete (header, key, 1); + snprintf (key, 10, "TFORM"OFF_T_FMT, zcol); + gfits_delete (header, key, 1); + + // create the output image + gfits_create_matrix (header, matrix); + + // 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]; + + // ztile[i] is the default (or max) tile size in the i-th dimension + // ntile[i] is the number of tiles in the i-th dimension + // otile[i] is the current output tile counter in the i-th dimension + } + + // this takes place in three major steps: + // 1) read the table data : zdef.format tells the size of the varlength heap element + // 2) uncompress the data : pixel size depends on compression method + // 3) distribute data to the tiles : output pixel size depends on header.bitpix + + // heapdata -> zdata -> odata -> idata + + // heapdata.pixsize : zdef.format + // zdata.pixsize : depends on compression + // odata.pixsize : depends on compression + // idata.pixsize : 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 + 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 final image (not needed) + // int tile_pixsize = abs(header[0].bitpix) / 8; + // fprintf (stderr, "raw_pixsize: %d, cmp_pixsize: %d, tile_pixsize: %d, raw_bitpix: %d\n", raw_pixsize, cmp_pixsize, tile_pixsize, raw_bitpix); + + int Nout_alloc = raw_pixsize*max_tile_size; + 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++) { + + // expected output size for this tile + Nout = raw_pixsize*gfits_tile_size (matrix, otile, ztile); + + zdata = gfits_varlength_column_pointer (ftable, &zdef, row, &Nzdata); + 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 : %d -> ", Nout); + } + + // 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 + 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, "%d\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"); + } + + if (!strcasecmp(cmptype, "GZIP_1")) { + // Nout is number of pixels + if (!gfits_byteswap_zdata (out, Nout * raw_pixsize, raw_pixsize)) ESCAPE; + } + + 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"); + } + + // copy the uncompressed pixels into their correct locations + if (!strcasecmp(cmptype, "GZIP_2") || !strcasecmp(cmptype, "NONE_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)) { + 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"); + } + + // 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; + } + } + } + + FREE (ztile); + if (optname != NULL) { + for (j = 0; j < Noptions; j++) { + FREE (optname[j]); + FREE (optvalue[j]); + } + FREE (optname); + FREE (optvalue); + } + FREE (out); + FREE (otile); + FREE (ntile); + return (TRUE); +} + +// bitpix is the input data size/type +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; + 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; + for (i = 1; i < matrix->Naxes; i++) { + Nline *= Ztile[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; + } + + // pixel offset in output array relative to tile start + int offset = 0; + + int directCopy = (zzero == 0.0) && (zscale == 1.0); + +# define SCALE_AND_DIST_INT_PRINT(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 = directCopy ? *RAWptr : *RAWptr * zscale + zzero; \ + fprintf (stderr, "dis: %08x : %08x\n", *TILEptr, *RAWptr); \ + } } } + + // 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_INT(OTYPE, OSIZE) { \ + OTYPE *TILEptr = (OTYPE *) &matrix->buffer[OSIZE*(offset + start)]; \ + for (j = 0; j < Ztile[0]; j++, TILEptr++, RAWptr++) { \ + if (*TILEptr == zblank) { \ + *TILEptr = oblank; \ + } else { \ + if (directCopy) { \ + *TILEptr = *RAWptr; \ + } 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(OTYPE, OSIZE) { \ + OTYPE *TILEptr = (OTYPE *) &matrix->buffer[OSIZE*(offset + start)]; \ + for (j = 0; j < Ztile[0]; j++, TILEptr++, RAWptr++) { \ + if (!isfinite(*TILEptr)) { \ + *TILEptr = oblank; \ + } else { \ + if (directCopy) { \ + *TILEptr = *RAWptr; \ + } else { \ + *TILEptr = *RAWptr * zscale + zzero; \ + } } } } + + // this macro sets up the outer switch and calls above macro with all output bitpix options +# define SETUP_RAWSIZE(ITYPE, ISIZE) { \ + ITYPE *RAWptr = (ITYPE *) &raw[i*ISIZE*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(); \ + } } + + // this macro sets up the outer switch and calls above macro with all output bitpix options +# define SETUP_RAWSIZE_PRINT(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 (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(); + } + + // 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]; + } + } + + free (counter); + free (Ztile); + return (TRUE); +} + +// 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*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: trunk/Ohana/src/libfits/matrix/F_uncompress_data.c =================================================================== --- trunk/Ohana/src/libfits/matrix/F_uncompress_data.c (revision 38255) +++ trunk/Ohana/src/libfits/matrix/F_uncompress_data.c (revision 38441) @@ -26,31 +26,42 @@ int gfits_uncompress (Bytef *dest, uLongf *destLen, const Bytef *source, uLong sourceLen); -int gfits_uncompress_data (char *zdata, int Nzdata, char *cmptype, char **optname, char **optvalue, int Nopt, char *outdata, int *Nout, int out_pixsize) { +# define ESCAPE(RET) { fprintf (stderr, "gzip error in %s @ %s:%d\n", __func__, __FILE__, __LINE__); return (RET); } + +int gfits_uncompress_data (char *zdata, int Nzdata, char *cmptype, char **optname, char **optvalue, int Nopt, char *outdata, int *Nout, int Nout_alloc, int out_pixsize) { int status; - static int Ninsum = 0; - static int Noutsum = 0; - if (!strcasecmp(cmptype, "GZIP_1")) { - unsigned long tNout = *Nout * out_pixsize; + // do not actually uncompress : this is used for testing + if (!strcasecmp(cmptype, "NONE") || !strcasecmp(cmptype, "NONE_1") || !strcasecmp(cmptype, "NONE_2")) { + memcpy (outdata, zdata, Nzdata); + myAssert (*Nout == Nzdata, "invalid size"); + *Nout = Nzdata / out_pixsize; + return (TRUE); + } + + if (!strcasecmp(cmptype, "GZIP_1") || !strcasecmp(cmptype, "GZIP_2")) { + // unsigned long tNout = *Nout * out_pixsize; + // input value of *Nout is number of BYTES + unsigned long tNout = Nout_alloc; // for GZIP data, I need to check for and remove the header on the first block: - gfits_gz_stripheader ((unsigned char *)zdata, &Nzdata); + // XXX maybe not anymore + // gfits_gz_stripheader ((unsigned char *)zdata, &Nzdata); // uncompress does not require us to know the expected number of pixel; it tells us the number + // XXX shouldn't we validate the result : we think we know the size status = gfits_uncompress ((Bytef *) outdata, &tNout, (Bytef *) zdata, Nzdata); - if (status != Z_OK) { - fprintf (stderr, "error in uncompress (GZIP)\n"); - return (FALSE); - } + if (status != Z_OK) ESCAPE(FALSE); + myAssert (*Nout == tNout, "uncompressed size mismatch"); *Nout = tNout / out_pixsize; + // output value of *Nout is number of PIXELS // the resulting uncompressed data is byteswapped - if (!gfits_byteswap_zdata (outdata, *Nout, out_pixsize)) return (FALSE); + // if (!gfits_byteswap_zdata (outdata, *Nout, out_pixsize)) return (FALSE); return (TRUE); } - if (!strcasecmp(cmptype, "RICE_1")) { + if (!strcasecmp(cmptype, "RICE_1") || !strcasecmp(cmptype, "RICE_ONE")) { int i, blocksize; // look for the BLOCKSIZE @@ -66,41 +77,36 @@ } - int Npix; - // Npix = *Nout * (out_pixsize / 4.0); - Npix = *Nout; - - Ninsum += Nzdata; - Noutsum += *Nout; - // fprintf (stderr, "%d comp bytes; %d uncomp 'pixels', totals: %d %d\n", Nzdata, Npix, Ninsum, Noutsum); + int status = FALSE; + int Npix = *Nout / out_pixsize; + int k; switch (out_pixsize) { 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); + status = fits_rdecomp ((unsigned char *) zdata, Nzdata, (unsigned int *) outdata, Npix, blocksize); + break; + case 2: + if (0) { + fprintf (stderr, "Nout: %d, Nrawpix: %d\n", Nzdata, Npix); + fprintf (stderr, "cmp out: "); + for (k = 0; k < Nzdata; k++) { fprintf (stderr, "0x%02hhx ", zdata[k]); } + fprintf (stderr, "\n"); } - return (TRUE); - - 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); - + status = fits_rdecomp_short ((unsigned char *) zdata, Nzdata, (unsigned short *) outdata, Npix, 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); - + status = fits_rdecomp_byte ((unsigned char *) zdata, Nzdata, (unsigned char *) outdata, Npix, blocksize); + break; default: fprintf (stderr, "invalid output pixel size %d\n", out_pixsize); - return (FALSE); + status = FALSE; + break; } - + if (status) { + *Nout = 0; + return FALSE; + } + *Nout = Npix; + return TRUE; } Index: trunk/Ohana/src/libfits/table/F_compress_T.c =================================================================== --- trunk/Ohana/src/libfits/table/F_compress_T.c (revision 38441) +++ trunk/Ohana/src/libfits/table/F_compress_T.c (revision 38441) @@ -0,0 +1,464 @@ +# include +# include +# include +# define VERBOSE_DUMP 0 + +// the user needs to specify the various compression options: +// ztilelen (may be 0) +// zcmptype + +// steps to convert an image to a compressed table: +// * determine the number of tiles (from ztile[] and image size) +// * construct an empty table with the right dimensions (1 column, Ntile rows) +// * loop over tiles +// * extract tile data to a buffer +// * compress the buffer data +// * insert in table heap +// * update header keywords + +// we do not try to support compression of an image with an associated table (not valid) + +# define ESCAPE { fprintf (stderr, "error in %s @ line %d\n", __func__, __LINE__); goto escape; } + +int gfits_collect_table_data (FTable *table, TableField *field, char *raw, int row_start, int Nrows); +int gfits_collect_table_gzp2 (FTable *table, TableField *field, char *raw, int row_start, int Nrows); + +static float timeSum1 = 0.0; +static float timeSum2 = 0.0; +static float timeSum2a = 0.0; +static float timeSum2b = 0.0; +static float timeSum2c = 0.0; +static float timeSum2d = 0.0; +static float timeSum2e = 0.0; +static float timeSum3 = 0.0; + +void gfits_compress_timing () { + fprintf (stderr, "cmp times: %f %f %f\n", timeSum1, timeSum2, timeSum3); + fprintf (stderr, "cmp times: %f %f %f %f %f\n", timeSum2a, timeSum2b, timeSum2c, timeSum2d, timeSum2e); + + timeSum1 = 0.0; + timeSum2 = 0.0; + timeSum2a = 0.0; + timeSum2b = 0.0; + timeSum2c = 0.0; + timeSum2d = 0.0; + timeSum2e = 0.0; + timeSum3 = 0.0; +} + +int gfits_compress_table (FTable *srctable, FTable *tgttable, int ztilelen, char *zcmptype) { + + int i, j; + + char keyword[81]; + + char *raw = NULL; + char *zdata = NULL; + + Header *srcheader = srctable->header; + Header *tgtheader = tgttable->header; + + // XXX EAM: + struct timeval startTimer, stopTimer; + float dtime; + gettimeofday (&startTimer, (void *) NULL); + + int Ntile, ztilelast; + if (!ztilelen) ztilelen = srcheader->Naxis[1]; + + if (!srcheader->Naxis[1]) { + Ntile = 0; + ztilelast = 0; + } else { + if (srcheader->Naxis[1] % ztilelen) { + Ntile = srcheader->Naxis[1] / ztilelen + 1; + ztilelast = srcheader->Naxis[1] % ztilelen; + } else { + Ntile = srcheader->Naxis[1] / ztilelen; + ztilelast = ztilelen; + } + } + + // EXTNAME is not mandatory, but I need to keep it + char extname[80]; + if (!gfits_scan (srcheader, "EXTNAME", "%s", 1, extname)) { + strcpy (extname, "COMPRESSED_TABLE"); + } + + // creates an empty (Naxes = 2, Naxis[i] = 0) table header with XTENSION = BINTABLE, EXTNAME = COMPRESSED_TABLE + if (!gfits_create_table_header (tgtheader, "BINTABLE", extname)) ESCAPE; + + int Nfields; + if (!gfits_scan (srcheader, "TFIELDS", "%d", 1, &Nfields)) ESCAPE; + + // XXX TIMER 1 + gettimeofday (&stopTimer, (void *) NULL); + dtime = DTIME (stopTimer, startTimer); + timeSum1 += dtime; + gettimeofday (&startTimer, (void *) NULL); + + TableField *fields = NULL; + ALLOCATE_ZERO (fields, TableField, Nfields); + for (i = 0; i < Nfields; i++) { + snprintf (keyword, 80, "TTYPE%d", i+1); + if (!gfits_scan (srcheader, keyword, "%s", 1, fields[i].ttype)) ESCAPE; + if (!gfits_scan_alt (srcheader, keyword, "%C", 1, fields[i].ttype_cmt)) ESCAPE; + + // TUNIT is not mandatory + snprintf (keyword, 80, "TUNIT%d", i+1); + if (!gfits_scan (srcheader, keyword, "%s", 1, fields[i].tunit)) { + fields[i].tunit[0] = 0; + } else { + if (!gfits_scan_alt (srcheader, keyword, "%C", 1, fields[i].tunit_cmt)) ESCAPE; + } + + snprintf (keyword, 80, "TFORM%d", i+1); + if (!gfits_scan (srcheader, keyword, "%s", 1, fields[i].tformat)) ESCAPE; + if (!gfits_scan_alt (srcheader, keyword, "%C", 1, fields[i].tformat_cmt)) ESCAPE; + + // for now we set all fields to the requested type. since I do not yet know the column data types I cannot yet assign automatic cmptypes + strcpy (fields[i].zctype, zcmptype); + + // by using "P" format, we are limited to 32bit pointers (2GB heap) + if (!gfits_define_bintable_column (tgtheader, "1QB(0)", fields[i].ttype, fields[i].ttype_cmt, fields[i].tunit, 1.0, 0.0)) ESCAPE; + } + + // allocates the default data array (tgttable->buffer) + if (!gfits_create_table (tgtheader, tgttable)) ESCAPE; + + // copy original header to output header (XXX this needs to be finished) + if (!gfits_copy_keywords_compress (srcheader, tgtheader)) ESCAPE; + + // this allocates the pointer data array for the full set of tiles (NOT the heap) + char *tmpbuffer = NULL; + ALLOCATE_ZERO (tmpbuffer, char, 16*Ntile); // need space for Ntile entries, each of width 16 bytes (2 long) + for (i = 0; i < Nfields; i++) { + if (!gfits_set_bintable_column (tgtheader, tgttable, fields[i].ttype, tmpbuffer, Ntile)) ESCAPE; + } + free (tmpbuffer); + + // add ZIMAGE from output header + if (!gfits_modify_alt (tgtheader, "ZTABLE", "%t", 1, TRUE)) ESCAPE; + + // define compression-specific keywords, update header as needed. + if (!gfits_modify (tgtheader, "ZTILELEN", "%d", 1, ztilelen)) ESCAPE; + + if (!gfits_modify (tgtheader, "ZNAXIS1", OFF_T_FMT, 1, srcheader->Naxis[0])) ESCAPE; + if (!gfits_modify (tgtheader, "ZNAXIS2", OFF_T_FMT, 1, srcheader->Naxis[1])) ESCAPE; + + off_t pcount; + if (!gfits_scan (srcheader, "PCOUNT", OFF_T_FMT, 1, &pcount)) ESCAPE; + if (!gfits_modify (tgtheader, "ZPCOUNT", OFF_T_FMT, 1, pcount)) ESCAPE; + + off_t theap; + if (gfits_scan (srcheader, "THEAP", OFF_T_FMT, 1, &theap)) { + if (!gfits_modify (tgtheader, "ZTHEAP", OFF_T_FMT, 1, theap)) ESCAPE; + } + + int max_width = 0; + int offset = 0; // bytes from first column of first field to the current field (accumulate to set) + + for (i = 0; i < Nfields; i++) { + snprintf (keyword, 81, "ZFORM%d", i+1); + if (!gfits_modify (tgtheader, keyword, "%s", 1, fields[i].tformat)) ESCAPE; + if (!gfits_varlength_column_define (tgttable, &fields[i].zdef, i+1)) ESCAPE; + if (!gfits_bintable_format (fields[i].tformat, fields[i].datatype, &fields[i].Nvalues, &fields[i].pixsize)) ESCAPE; // + fields[i].rowsize = fields[i].Nvalues*fields[i].pixsize; + max_width = MAX(max_width, fields[i].rowsize); + + if (!strcasecmp (fields[i].zctype, "AUTO")) { + if (!strcmp (fields[i].datatype, "short") || + !strcmp (fields[i].datatype, "float") || + !strcmp (fields[i].datatype, "double")|| + !strcmp (fields[i].datatype, "int64_t")) { + strcpy (fields[i].zctype, "GZIP_2"); + goto got_cmptype; + } + if (!strcmp (fields[i].datatype, "int")) { + strcpy (fields[i].zctype, "RICE_1"); + goto got_cmptype; + } + if (!strcmp (fields[i].datatype, "byte") || + !strcmp (fields[i].datatype, "char")) { + strcpy (fields[i].zctype, "GZIP_1"); + goto got_cmptype; + } + } + got_cmptype: + + + // OVERRIDE: RICE can only be used on integer fields + if (!strcasecmp (fields[i].zctype, "RICE_1") || !strcasecmp (fields[i].zctype, "RICE_ONE")) { + if (!strcmp (fields[i].datatype, "float") || !strcmp (fields[i].datatype, "double") || !strcmp (fields[i].datatype, "int64_t")) { + strcpy (fields[i].zctype, "GZIP_2"); + } + } + + // OVERRIDE: GZIP_2 invalid for B (use GZIP_1) + if (!strcasecmp (fields[i].zctype, "GZIP_2") && !strcmp (fields[i].datatype, "byte")) { + strcpy (fields[i].zctype, "GZIP_1"); + } + + // compression type per column (XXX for now we are just using zcmptype, as set above) + snprintf (keyword, 81, "ZCTYP%d", i+1); + if (!gfits_modify (tgtheader, keyword, "%s", 1, fields[i].zctype)) ESCAPE; + + fields[i].offset = offset; + offset += fields[i].rowsize; + } + + // allocate the intermediate storage buffers + int Nzdata_alloc = max_width*ztilelen + 100; + ALLOCATE (raw, char, max_width*ztilelen); + ALLOCATE (zdata, char, Nzdata_alloc); + + gfits_dump_raw_table (srctable, "cmp"); + + // XXX TIMER 2 + gettimeofday (&stopTimer, (void *) NULL); + dtime = DTIME (stopTimer, startTimer); + timeSum2 += dtime; + gettimeofday (&startTimer, (void *) NULL); + + // compress the data : copy into a tile, compress the tile, then add to the output table + // each tile -> 1 row of the output table + for (i = 0; i < Ntile; i++) { + + for (j = 0; j < Nfields; j++) { + + gettimeofday (&startTimer, (void *) NULL); + + int Nrows = (i == Ntile - 1) ? ztilelast : ztilelen; + int row_start = i*ztilelen; + // ^- first row for this tile & field + // NOTE: assumes each tile has the same length, ex the last (true for now) + + // copy the raw pixels from their native matrix locations to the temporary output buffer + if (!strcasecmp(fields[j].zctype, "GZIP_2") || !strcasecmp(fields[j].zctype, "NONE_2")) { + if (!gfits_collect_table_gzp2 (srctable, &fields[j], raw, row_start, Nrows)) ESCAPE; + } else { + if (!gfits_collect_table_data (srctable, &fields[j], raw, row_start, Nrows)) ESCAPE; + } + + // XXX TIMER 2a + gettimeofday (&stopTimer, (void *) NULL); + dtime = DTIME (stopTimer, startTimer); + timeSum2a += dtime; + gettimeofday (&startTimer, (void *) NULL); + + if (VERBOSE_DUMP) { + int k; + fprintf (stderr, "c1: "); + for (k = 0; k < Nrows*fields[j].Nvalues*fields[j].pixsize; k++) { + fprintf (stderr, "%02hhx", raw[k]); + if (k % 2) fprintf (stderr, " "); + // if (k % 32 == 31) fprintf (stderr, "\n"); + } + fprintf (stderr, "\n"); + } + + // optname, optvalue = NULL, Noptions = 0 + + int Nraw = Nrows*fields[j].Nvalues; // number of pixels + if (!strcasecmp(fields[j].zctype, "GZIP_1")) { + if (!gfits_byteswap_zdata (raw, Nraw * fields[j].pixsize, fields[j].pixsize)) ESCAPE; + } + + // XXX TIMER 2b + gettimeofday (&stopTimer, (void *) NULL); + dtime = DTIME (stopTimer, startTimer); + timeSum2b += dtime; + gettimeofday (&startTimer, (void *) NULL); + + if (VERBOSE_DUMP) { + int k; + fprintf (stderr, "c2: "); + for (k = 0; k < Nrows*fields[j].Nvalues*fields[j].pixsize; k++) { + fprintf (stderr, "%02hhx", raw[k]); + if (k % 2) fprintf (stderr, " "); + // if (k % 32 == 31) fprintf (stderr, "\n"); + } + fprintf (stderr, "\n"); + } + + 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, 0, 0)) ESCAPE; + + // XXX TIMER 2c + gettimeofday (&stopTimer, (void *) NULL); + dtime = DTIME (stopTimer, startTimer); + timeSum2c += dtime; + gettimeofday (&startTimer, (void *) NULL); + + if (VERBOSE_DUMP) { + int k; + fprintf (stderr, "c3: "); + for (k = 0; k < Nzdata; k++) { + fprintf (stderr, "%02hhx", zdata[k]); + if (k % 2) fprintf (stderr, " "); + // if (k % 32 == 31) fprintf (stderr, "\n"); + } + fprintf (stderr, "\n"); + } + + if (strcasecmp(fields[j].zctype, "NONE_2") && // NONE and NONE_1 not swapped? + strcasecmp(fields[j].zctype, "GZIP_1") && + strcasecmp(fields[j].zctype, "GZIP_2") && + strcasecmp(fields[j].zctype, "RICE_1") && + strcasecmp(fields[j].zctype, "RICE_ONE")) { + if (!gfits_byteswap_zdata (zdata, Nzdata, fields[j].pixsize)) ESCAPE; + } + + // XXX TIMER 2d + gettimeofday (&stopTimer, (void *) NULL); + dtime = DTIME (stopTimer, startTimer); + timeSum2d += dtime; + gettimeofday (&startTimer, (void *) NULL); + + if (VERBOSE_DUMP) { + int k; + fprintf (stderr, "c4: "); + for (k = 0; k < Nzdata; k++) { + fprintf (stderr, "%02hhx", zdata[k]); + if (k % 2) fprintf (stderr, " "); + // if (k % 32 == 31) fprintf (stderr, "\n"); + } + fprintf (stderr, "\n"); + } + + if (!gfits_varlength_column_add_data (tgttable, zdata, Nzdata, i, &fields[j].zdef)) ESCAPE; + // XXX TIMER 2e + gettimeofday (&stopTimer, (void *) NULL); + dtime = DTIME (stopTimer, startTimer); + timeSum2e += dtime; + gettimeofday (&startTimer, (void *) NULL); + + } + } + + for (i = 0; i < Nfields; i++) { + if (!gfits_varlength_column_finish (tgttable, &fields[i].zdef)) ESCAPE; + } + + // XXX TIMER 3 + gettimeofday (&stopTimer, (void *) NULL); + dtime = DTIME (stopTimer, startTimer); + timeSum3 += dtime; + gettimeofday (&startTimer, (void *) NULL); + + gfits_dump_cmp_table (tgttable, "cmp"); + + free (raw); + free (zdata); + return (TRUE); + +escape: + FREE (raw); + FREE (zdata); + return FALSE; +} + +# define VERBOSE 0 + +// 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_table_data (FTable *table, TableField *field, char *raw, int row_start, int Nrows) { + + off_t i; + + // we are copying NN rows into the column which starts at XX and has MM bytes per row + + off_t Nx = table->header->Naxis[0]; + + char *tblbuffer = &table->buffer[Nx*row_start + field->offset]; + + int rowsize = field->rowsize; + + if (VERBOSE) fprintf (stderr, "distribute: "); + for (i = 0; i < Nrows; i++, tblbuffer += Nx) { + memcpy (&raw[i*rowsize], tblbuffer, rowsize); +# if (VERBOSE) + int j; for (j = 0; j < field->rowsize; j++) fprintf (stderr, "0x%02hhx ", table->buffer[Nx*row + field->offset + j]); +# endif + } + if (VERBOSE) fprintf (stderr, "\n"); + + return (TRUE); +} + +int gfits_collect_table_data_alt (FTable *table, TableField *field, char *raw, int row_start, int Nrows) { + + off_t i; + + // we are copying NN rows of the column which starts at XX and has MM bytes per row + + off_t Nx = table->header->Naxis[0]; + + if (VERBOSE) fprintf (stderr, "collect: "); + for (i = 0; i < Nrows; i++) { + int row = row_start + i; + memcpy (&raw[i*field->rowsize], &table->buffer[Nx*row + field->offset], field->rowsize); +# if (VERBOSE) + int j; for (j = 0; j < field->rowsize; j++) fprintf (stderr, "0x%02hhx ", table->buffer[Nx*row + field->offset + j]); +# endif + } + if (VERBOSE) fprintf (stderr, "\n"); + + return (TRUE); +} + +int gfits_collect_table_gzp2 (FTable *table, TableField *field, char *raw, int row_start, int Nrows) { + + off_t i, j, k; + + // we are copying NN rows into the column which starts at XX and has MM bytes per row + + off_t Nx = table->header->Naxis[0]; + int Nvalues = field->Nvalues; + int pixsize = field->pixsize; + int offset = field->offset; + + char *rawptr = raw; + + for (k = 0; k < field->pixsize; k++) { +# ifdef BYTE_SWAP + char *tblptr_start = &table->buffer[Nx*row_start + offset + (pixsize - k - 1)]; +# else + char *tblptr_start = &table->buffer[Nx*row_start + offset + k]; +# endif + for (i = 0; i < Nrows; i++, tblptr_start += Nx) { + char *tblptr = tblptr_start; + for (j = 0; j < Nvalues; j++, tblptr += pixsize, rawptr++) { + *rawptr = *tblptr; + } + } + } + return (TRUE); +} + +int gfits_collect_table_gzp2_alt (FTable *table, TableField *field, char *raw, int row_start, int Nrows) { + + off_t i, j, k; + + // we are copying NN rows into the column which starts at XX and has MM bytes per row + + off_t Nx = table->header->Naxis[0]; + + for (k = 0; k < field->pixsize; k++) { + for (i = 0; i < Nrows; i++) { + int row = row_start + i; + char *rawptr = &raw[i*field->Nvalues + k*Nrows*field->Nvalues]; +# ifdef BYTE_SWAP + char *tblptr = &table->buffer[Nx*row + field->offset + (field->pixsize - k - 1)]; +# else + char *tblptr = &table->buffer[Nx*row + field->offset + k]; +# endif + for (j = 0; j < field->Nvalues; j++, tblptr += field->pixsize, rawptr++) { + *rawptr = *tblptr; + myAssert (rawptr - raw < Nrows*field->Nvalues*field->pixsize, "oops"); + myAssert (tblptr - table->buffer < table->header->Naxis[0]*table->header->Naxis[1], "oops"); + } + } + } + return (TRUE); +} Index: trunk/Ohana/src/libfits/table/F_copy_T.c =================================================================== --- trunk/Ohana/src/libfits/table/F_copy_T.c (revision 38255) +++ trunk/Ohana/src/libfits/table/F_copy_T.c (revision 38441) @@ -6,9 +6,25 @@ /* find buffer size */ - out[0].validsize = in[0].validsize; - out[0].datasize = in[0].datasize; + out[0].validsize = in[0].validsize; + out[0].datasize = in[0].datasize; + out[0].heap_start = in[0].heap_start; ALLOCATE (out[0].buffer, char, out[0].datasize); memcpy (out[0].buffer, in[0].buffer, out[0].datasize); + return (TRUE); +} + +/*********************** copy ftable pointers (does not copy memory) ***********************************/ +int gfits_copy_ftable_ptr (FTable *in, FTable *out) { + + if (!in) return FALSE; + if (!out) return FALSE; + + /* find buffer size */ + out[0].validsize = in[0].validsize; + out[0].datasize = in[0].datasize; + out[0].heap_start = in[0].heap_start; + out[0].header = in[0].header; + out[0].buffer = in[0].buffer; return (TRUE); } @@ -28,7 +44,8 @@ /* find buffer size */ - out[0].datasize = in[0].datasize; - out[0].pad = in[0].pad; - out[0].Nrow = in[0].Nrow; + out[0].datasize = in[0].datasize; + // out[0].heap_start = in[0].heap_start; + out[0].pad = in[0].pad; + out[0].Nrow = in[0].Nrow; off_t Nrows = out[0].Nrow; Index: trunk/Ohana/src/libfits/table/F_create_T.c =================================================================== --- trunk/Ohana/src/libfits/table/F_create_T.c (revision 38255) +++ trunk/Ohana/src/libfits/table/F_create_T.c (revision 38441) @@ -4,8 +4,9 @@ /*********************** fits init table *******************************/ int gfits_init_table (FTable *table) { - table[0].header = NULL; - table[0].buffer = NULL; - table[0].datasize = 0; - table[0].validsize = 0; + table[0].header = NULL; + table[0].buffer = NULL; + table[0].datasize = 0; + table[0].heap_start = 0; + table[0].validsize = 0; return (TRUE); @@ -42,6 +43,21 @@ } table[0].datasize = Nbytes; + table[0].heap_start = gfits_heap_start (header); return (TRUE); +} + +/*********************** fits create table *******************************/ +int gfits_set_table_rows (Header *header, FTable *table, off_t Nrows) { + + header[0].Naxis[1] = Nrows; + gfits_modify (header, "NAXIS2", OFF_T_FMT, 1, Nrows); + + off_t Nbytes = gfits_data_size (header); + REALLOCATE (table[0].buffer, char, Nbytes); + bzero (table[0].buffer, Nbytes); + table[0].datasize = Nbytes; + table[0].heap_start = gfits_heap_start (header); + return (TRUE); } @@ -60,4 +76,5 @@ Ny = table[0].header[0].Naxis[0]; table[0].datasize = gfits_data_size (header); + // table[0].heap_start = gfits_heap_start (header); table[0].pad = table[0].datasize - Nx*Ny; Index: trunk/Ohana/src/libfits/table/F_define_column.c =================================================================== --- trunk/Ohana/src/libfits/table/F_define_column.c (revision 38255) +++ trunk/Ohana/src/libfits/table/F_define_column.c (revision 38441) @@ -66,7 +66,11 @@ sprintf (field, "TTYPE%d", Nfields); gfits_modify (header, field, "%s", 1, label); - gfits_modify_alt (header, field, "%C", 1, comment); - sprintf (field, "TUNIT%d", Nfields); - gfits_modify (header, field, "%s", 1, unit); + if (comment) { + gfits_modify_alt (header, field, "%C", 1, comment); + } + if (unit) { + sprintf (field, "TUNIT%d", Nfields); + gfits_modify (header, field, "%s", 1, unit); + } sprintf (field, "TFORM%d", Nfields); gfits_modify (header, field, "%s", 1, format); Index: trunk/Ohana/src/libfits/table/F_get_column.c =================================================================== --- trunk/Ohana/src/libfits/table/F_get_column.c (revision 38255) +++ trunk/Ohana/src/libfits/table/F_get_column.c (revision 38441) @@ -17,4 +17,9 @@ void *gfits_get_bintable_column_data (Header *header, FTable *table, char *label, char *type, off_t *Nrow, int *Ncol) { + void *data = gfits_get_bintable_column_data_raw (header, table, label, type, Nrow, Ncol, FALSE); + return data; +} + +void *gfits_get_bintable_column_data_raw (Header *header, FTable *table, char *label, char *type, off_t *Nrow, int *Ncol, char nativeOrder) { off_t Nx, Ny; @@ -71,15 +76,19 @@ } - /* convert data in-situ with correct type, byte swap and Bzero/Bscale */ + // NOTE: we have already copied the data to 'array', so the blocks below + // only need to swap if this is a direct copy Pin = array; Pout = array; + int directCopy = (Bzero == 0.0) && (Bscale == 1.0); + + /* convert data in-situ with correct type, byte swap and Bzero/Bscale */ if (!strcmp (type, "char")) { for (i = 0; i < Nval*Ny; i++, Pin+=Nbytes, Pout+=Nbytes) { - *(char *)Pout = *(char *)Pin*Bscale + Bzero; + if (!directCopy) { *(char *)Pout = *(char *)Pin*Bscale + Bzero; } } } if (!strcmp (type, "byte")) { for (i = 0; i < Nval*Ny; i++, Pin+=Nbytes, Pout+=Nbytes) { - *(char *)Pout = *(char *)Pin*Bscale + Bzero; + if (!directCopy) { *(char *)Pout = *(char *)Pin*Bscale + Bzero; } } } @@ -87,7 +96,7 @@ for (i = 0; i < Nval*Ny; i++, Pin+=Nbytes, Pout+=Nbytes) { # ifdef BYTE_SWAP - SWAP_BYTE; -# endif - *(short *)Pout = *(short *)Pin*Bscale + Bzero; + if (!nativeOrder) { SWAP_BYTE; } +# endif + if (!directCopy) { *(short *)Pout = *(short *)Pin*Bscale + Bzero; } } } @@ -95,24 +104,15 @@ for (i = 0; i < Nval*Ny; i++, Pin+=Nbytes, Pout+=Nbytes) { # ifdef BYTE_SWAP - SWAP_WORD; -# endif - *(int *)Pout = *(int *)Pin*Bscale + Bzero; + if (!nativeOrder) { SWAP_WORD; } +# endif + if (!directCopy) { *(int *)Pout = *(int *)Pin*Bscale + Bzero; } } } if (!strcmp (type, "int64_t")) { - if ((Bzero == 0.0) && (Bscale == 1.0)) { - for (i = 0; i < Nval*Ny; i++, Pin+=Nbytes, Pout+=Nbytes) { -# ifdef BYTE_SWAP - SWAP_DBLE; -# endif - *(int64_t *)Pout = *(int64_t *)Pin; - } - } else { - for (i = 0; i < Nval*Ny; i++, Pin+=Nbytes, Pout+=Nbytes) { -# ifdef BYTE_SWAP - SWAP_DBLE; -# endif - *(int64_t *)Pout = *(int64_t *)Pin*Bscale + Bzero; - } + for (i = 0; i < Nval*Ny; i++, Pin+=Nbytes, Pout+=Nbytes) { +# ifdef BYTE_SWAP + if (!nativeOrder) { SWAP_DBLE; } +# endif + if (!directCopy) { *(int64_t *)Pout = *(int64_t *)Pin*Bscale + Bzero; } } } @@ -120,7 +120,7 @@ for (i = 0; i < Nval*Ny; i++, Pin+=Nbytes, Pout+=Nbytes) { # ifdef BYTE_SWAP - SWAP_WORD; -# endif - *(float *)Pout = *(float *)Pin*Bscale + Bzero; + if (!nativeOrder) { SWAP_WORD; } +# endif + if (!directCopy) { *(float *)Pout = *(float *)Pin*Bscale + Bzero; } } } @@ -128,9 +128,11 @@ for (i = 0; i < Nval*Ny; i++, Pin+=Nbytes, Pout+=Nbytes) { # ifdef BYTE_SWAP - SWAP_DBLE; -# endif - *(double *)Pout = *(double *)Pin*Bscale + Bzero; - } - } + if (!nativeOrder) { SWAP_DBLE; } +# endif + if (!directCopy) { *(double *)Pout = *(double *)Pin*Bscale + Bzero; } + } + } + + // check that we supplied a valid type *Ncol = Nval; @@ -174,4 +176,18 @@ if (!gfits_get_bintable_column_type_by_N (header, N, type, Nval)) return FALSE; return (TRUE); +} + +/***********************/ +int gfits_get_bintable_column_raw (Header *header, FTable *table, char *label, void **data, char nativeOrder) { + + char type[16]; + off_t Nrow; + int Ncol; + + char *array = gfits_get_bintable_column_data_raw (header, table, label, type, &Nrow, &Ncol, nativeOrder); + if (array == NULL) return (FALSE); + + *data = array; + return TRUE; } Index: trunk/Ohana/src/libfits/table/F_read_T.c =================================================================== --- trunk/Ohana/src/libfits/table/F_read_T.c (revision 38255) +++ trunk/Ohana/src/libfits/table/F_read_T.c (revision 38441) @@ -45,19 +45,17 @@ if (Nread < gfits_data_min_size (table[0].header)) { fprintf (stderr, "error: fits read error in %s, read "OFF_T_FMT", need "OFF_T_FMT"\n", __func__, Nread, gfits_data_min_size (table[0].header)); - if (padIfShort) { - memset (&table[0].buffer[Nread], 0, Nbytes - Nread); - fprintf (stderr, "warning: file missing data, padding with zeros: USE AT YOUR OWN RISK!\n"); - table[0].validsize = Nread; - table[0].datasize = Nbytes; - return (TRUE); - } else { - gfits_free_table (table); + if (!padIfShort) { + gfits_free_table (table); return (FALSE); } - } - fprintf (stderr, "warning: file missing pad\n"); + memset (&table[0].buffer[Nread], 0, Nbytes - Nread); + fprintf (stderr, "warning: file missing data, padding with zeros: USE AT YOUR OWN RISK!\n"); + } else { + fprintf (stderr, "warning: file missing pad\n"); + } } table[0].validsize = Nread; table[0].datasize = Nbytes; + table[0].heap_start = gfits_heap_start (table->header); return (TRUE); } @@ -93,13 +91,10 @@ } fprintf (stderr, "error: fits read error in %s, read "OFF_T_FMT", need "OFF_T_FMT"\n", __func__, Nread, Nbytes); - if (padIfShort) { - memset (&table[0].buffer[Nread], 0, Nbytes - Nread); - fprintf (stderr, "warning: file missing data, padding with zeros: USE AT YOUR OWN RISK!\n"); - table[0].validsize = Nread; - table[0].datasize = Nbytes; - } else { + if (!padIfShort) { free (table[0].buffer); return (FALSE); } + memset (&table[0].buffer[Nread], 0, Nbytes - Nread); + fprintf (stderr, "warning: file missing data, padding with zeros: USE AT YOUR OWN RISK!\n"); } @@ -107,5 +102,7 @@ table[0].header[0].Naxis[1] = Nrows; gfits_modify (table[0].header, "NAXIS2", OFF_T_FMT, 1, Nrows); - table[0].datasize = gfits_data_size (table[0].header); + table[0].validsize = Nread; + table[0].datasize = Nbytes; + table[0].heap_start = gfits_heap_start (table->header); return (TRUE); @@ -122,4 +119,5 @@ Ny = table[0].header[0].Naxis[1]; table[0].datasize = gfits_data_size (table[0].header); + // table[0].heap_start = gfits_heap_start (table[0].header); table[0].pad = table[0].datasize - Nx*Ny; @@ -206,4 +204,5 @@ for (i = 0; i < Nrow; i++) table[0].row[i] = row[i]; table[0].datasize = gfits_data_size (table[0].header); + // table[0].heap_start = gfits_heap_start (table[0].header); table[0].pad = table[0].datasize - Nx*Ny; return (TRUE); Index: trunk/Ohana/src/libfits/table/F_set_column.c =================================================================== --- trunk/Ohana/src/libfits/table/F_set_column.c (revision 38255) +++ trunk/Ohana/src/libfits/table/F_set_column.c (revision 38441) @@ -1,4 +1,8 @@ # include # include + +# define SWAP_NONE + +# ifdef BYTE_SWAP # define SWAP_BYTE { \ char tmp; \ @@ -14,9 +18,14 @@ tmp = Pout[2]; Pout[2] = Pout[5]; Pout[5] = tmp; \ tmp = Pout[3]; Pout[3] = Pout[4]; Pout[4] = tmp; } +# else +# define SWAP_BYTE +# define SWAP_WORD +# define SWAP_DBLE +# endif /***********************/ int gfits_set_bintable_column (Header *header, FTable *table, char *label, void *data, off_t Nrow) { - off_t Nx, Ny, nbytes; + off_t Nx, Ny; int i, N, Nfields; int Nval, Nbytes, Nstart, Nv, Nb; @@ -55,13 +64,9 @@ gfits_scan (header, "NAXIS1", OFF_T_FMT, 1, &Nx); gfits_scan (header, "NAXIS2", OFF_T_FMT, 1, &Ny); + + // if no rows have yet been assigned, we need to allocate the full data buffer if (Ny == 0) { Ny = Nrow; - header[0].Naxis[1] = Ny; - gfits_modify (header, "NAXIS2", OFF_T_FMT, 1, Ny); - - nbytes = gfits_data_size (header); - REALLOCATE (table[0].buffer, char, MAX (nbytes, 1)); - bzero (table[0].buffer, nbytes); - table[0].datasize = nbytes; + gfits_set_table_rows (header, table, Ny); } if (Ny != Nrow) return (FALSE); @@ -143,4 +148,22 @@ } + if (!strcmp (type, "var")) { + for (i = 0; i < Nval*Nrow; i++, Pin+=Nbytes, Pout+=Nbytes) { + *(int *)Pout = *(int *)Pin; +# ifdef BYTE_SWAP + SWAP_WORD; +# endif + } + } + + if (!strcmp (type, "var64")) { + for (i = 0; i < Nval*Nrow; i++, Pin+=Nbytes, Pout+=Nbytes) { + *(int *)Pout = *(int *)Pin; +# ifdef BYTE_SWAP + SWAP_DBLE; +# endif + } + } + /* check array space */ if (Nx*Ny < Nx*(Nrow - 1) + Nstart + Nval*Nbytes) { @@ -161,7 +184,8 @@ /***********************/ -int gfits_set_bintable_column_reformat (Header *header, FTable *table, char *label, char *intype, void *data, off_t Nrow) { - - off_t Nx, Ny, nbytes; +// convert the input data array (of the specified intype) to the desired table data type. swap unless nativeOrder is requested +int gfits_set_bintable_column_reformat (Header *header, FTable *table, char *label, char *intype, void *data, off_t Nrow, char nativeOrder) { + + off_t Nx, Ny; int i, N, Nfields; int Nval, NbytesOut, Nstart, Nv, Nb; @@ -202,11 +226,5 @@ if (Ny == 0) { Ny = Nrow; - header[0].Naxis[1] = Ny; - gfits_modify (header, "NAXIS2", OFF_T_FMT, 1, Ny); - - nbytes = gfits_data_size (header); - REALLOCATE (table[0].buffer, char, MAX (nbytes, 1)); - bzero (table[0].buffer, nbytes); - table[0].datasize = nbytes; + gfits_set_table_rows (header, table, Ny); } if (Ny != Nrow) return (FALSE); @@ -230,12 +248,141 @@ // if (!strcmp (intype, #ITYPE)) { + int directCopy = (Bzero == 0.0) && (Bscale == 1.0); + +# define SET_VALUES(OUTNAME, OUTTYPE, INNAME, INTYPE, SWAP_OP, NBYTES_IN) \ + if (!strcmp (outtype, OUTNAME) && !strcmp (intype, INNAME)) { \ + int NbytesIn = NBYTES_IN; \ + for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) { \ + if (directCopy) { \ + *(OUTTYPE *)Pout = *(INTYPE *)Pin; \ + } else { \ + *(OUTTYPE *)Pout = (*(INTYPE *)Pin - Bzero) / Bscale; \ + } \ + if (!nativeOrder) { SWAP_OP; }}} + + SET_VALUES("char", char, "char", char, SWAP_NONE, 1); + SET_VALUES("byte", char, "char", char, SWAP_NONE, 1); + SET_VALUES("short", short, "char", char, SWAP_BYTE, 1); + SET_VALUES("int", int, "char", char, SWAP_WORD, 1); + SET_VALUES("int64_t", int64_t, "char", char, SWAP_DBLE, 1); + SET_VALUES("float", float, "char", char, SWAP_WORD, 1); + SET_VALUES("double", double, "char", char, SWAP_DBLE, 1); + + SET_VALUES("char", char, "byte", char, SWAP_NONE, 1); + SET_VALUES("byte", char, "byte", char, SWAP_NONE, 1); + SET_VALUES("short", short, "byte", char, SWAP_BYTE, 1); + SET_VALUES("int", int, "byte", char, SWAP_WORD, 1); + SET_VALUES("int64_t", int64_t, "byte", char, SWAP_DBLE, 1); + SET_VALUES("float", float, "byte", char, SWAP_WORD, 1); + SET_VALUES("double", double, "byte", char, SWAP_DBLE, 1); + + SET_VALUES("char", char, "short", short, SWAP_NONE, 2); + SET_VALUES("byte", char, "short", short, SWAP_NONE, 2); + SET_VALUES("short", short, "short", short, SWAP_BYTE, 2); + SET_VALUES("int", int, "short", short, SWAP_WORD, 2); + SET_VALUES("int64_t", int64_t, "short", short, SWAP_DBLE, 2); + SET_VALUES("float", float, "short", short, SWAP_WORD, 2); + SET_VALUES("double", double, "short", short, SWAP_DBLE, 2); + + SET_VALUES("char", char, "int", int, SWAP_NONE, 4); + SET_VALUES("byte", char, "int", int, SWAP_NONE, 4); + SET_VALUES("short", short, "int", int, SWAP_BYTE, 4); + SET_VALUES("int", int, "int", int, SWAP_WORD, 4); + SET_VALUES("int64_t", int64_t, "int", int, SWAP_DBLE, 4); + SET_VALUES("float", float, "int", int, SWAP_WORD, 4); + SET_VALUES("double", double, "int", int, SWAP_DBLE, 4); + + SET_VALUES("char", char, "int64_t", int64_t, SWAP_NONE, 8); + SET_VALUES("byte", char, "int64_t", int64_t, SWAP_NONE, 8); + SET_VALUES("short", short, "int64_t", int64_t, SWAP_BYTE, 8); + SET_VALUES("int", int, "int64_t", int64_t, SWAP_WORD, 8); + SET_VALUES("int64_t", int64_t, "int64_t", int64_t, SWAP_DBLE, 8); + SET_VALUES("float", float, "int64_t", int64_t, SWAP_WORD, 8); + SET_VALUES("double", double, "int64_t", int64_t, SWAP_DBLE, 8); + + SET_VALUES("char", char, "float", float, SWAP_NONE, 4); + SET_VALUES("byte", char, "float", float, SWAP_NONE, 4); + SET_VALUES("short", short, "float", float, SWAP_BYTE, 4); + SET_VALUES("int", int, "float", float, SWAP_WORD, 4); + SET_VALUES("int64_t", int64_t, "float", float, SWAP_DBLE, 4); + SET_VALUES("float", float, "float", float, SWAP_WORD, 4); + SET_VALUES("double", double, "float", float, SWAP_DBLE, 4); + + SET_VALUES("char", char, "double", double, SWAP_NONE, 8); + SET_VALUES("byte", char, "double", double, SWAP_NONE, 8); + SET_VALUES("short", short, "double", double, SWAP_BYTE, 8); + SET_VALUES("int", int, "double", double, SWAP_WORD, 8); + SET_VALUES("int64_t", int64_t, "double", double, SWAP_DBLE, 8); + SET_VALUES("float", float, "double", double, SWAP_WORD, 8); + SET_VALUES("double", double, "double", double, SWAP_DBLE, 8); + +# if (0) /** input == char **/ if (!strcmp (outtype, "char") && !strcmp (intype, "char")) { int NbytesIn = 1; for (i = 0; i < Nval*Nrow; i++, Pin+=NbytesIn, Pout+=NbytesOut) { - *(char *)Pout = (*(char *)Pin - Bzero) / Bscale; + if (directCopy) { *(char *)Pout = *(char *)Pin; } else { + *(char *)Pout = (*(char *)Pin - Bzero) / Bscale; + } } } if (!strcmp (outtype, "byte") && !strcmp (intype, "char")) { + int NbytesIn = 1; + for (i = 0; i < Nval*Nrow; i++, Pin+=NbytesIn, Pout+=NbytesOut) { + if (directCopy) { *(char *)Pout = *(char *)Pin; } else { + *(char *)Pout = (*(char *)Pin - Bzero) / Bscale; + } + } + } + if (!strcmp (outtype, "short") && !strcmp (intype, "char")) { + int NbytesIn = 1; + for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) { + if (directCopy) { *(short *)Pout = *(char *)Pin; } else { + *(short *)Pout = (*(char *)Pin - Bzero) / Bscale; + } +# ifdef BYTE_SWAP + if (!nativeOrder) { SWAP_BYTE; } +# endif + } + } + if (!strcmp (outtype, "int") && !strcmp (intype, "char")) { + int NbytesIn = 1; + for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) { + *(int *)Pout = (*(char *)Pin - Bzero) / Bscale; +# ifdef BYTE_SWAP + if (!nativeOrder) { SWAP_WORD; } +# endif + } + } + if (!strcmp (outtype, "int64_t") && !strcmp (intype, "char")) { + int NbytesIn = 1; + for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) { + *(int64_t *)Pout = (*(char *)Pin - Bzero) / Bscale; +# ifdef BYTE_SWAP + if (!nativeOrder) { SWAP_DBLE; } +# endif + } + } + if (!strcmp (outtype, "float") && !strcmp (intype, "char")) { + int NbytesIn = 1; + for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) { + *(float *)Pout = (*(char *)Pin - Bzero) / Bscale; +# ifdef BYTE_SWAP + if (!nativeOrder) { SWAP_WORD; } +# endif + } + } + if (!strcmp (outtype, "double") && !strcmp (intype, "char")) { + int NbytesIn = 1; + for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) { + *(double *)Pout = (*(char *)Pin - Bzero) / Bscale; +# ifdef BYTE_SWAP + if (!nativeOrder) { SWAP_DBLE; } +# endif + } + } + + /** input == byte **/ + if (!strcmp (outtype, "char") && !strcmp (intype, "byte")) { int NbytesIn = 1; for (i = 0; i < Nval*Nrow; i++, Pin+=NbytesIn, Pout+=NbytesOut) { @@ -243,52 +390,5 @@ } } - if (!strcmp (outtype, "short") && !strcmp (intype, "char")) { - int NbytesIn = 1; - for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) { - *(short *)Pout = (*(char *)Pin - Bzero) / Bscale; -# ifdef BYTE_SWAP - SWAP_BYTE; -# endif - } - } - if (!strcmp (outtype, "int") && !strcmp (intype, "char")) { - int NbytesIn = 1; - for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) { - *(int *)Pout = (*(char *)Pin - Bzero) / Bscale; -# ifdef BYTE_SWAP - SWAP_WORD; -# endif - } - } - if (!strcmp (outtype, "int64_t") && !strcmp (intype, "char")) { - int NbytesIn = 1; - for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) { - *(int64_t *)Pout = (*(char *)Pin - Bzero) / Bscale; -# ifdef BYTE_SWAP - SWAP_DBLE; -# endif - } - } - if (!strcmp (outtype, "float") && !strcmp (intype, "char")) { - int NbytesIn = 1; - for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) { - *(float *)Pout = (*(char *)Pin - Bzero) / Bscale; -# ifdef BYTE_SWAP - SWAP_WORD; -# endif - } - } - if (!strcmp (outtype, "double") && !strcmp (intype, "char")) { - int NbytesIn = 1; - for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) { - *(double *)Pout = (*(char *)Pin - Bzero) / Bscale; -# ifdef BYTE_SWAP - SWAP_DBLE; -# endif - } - } - - /** input == byte **/ - if (!strcmp (outtype, "char") && !strcmp (intype, "byte")) { + if (!strcmp (outtype, "byte") && !strcmp (intype, "byte")) { int NbytesIn = 1; for (i = 0; i < Nval*Nrow; i++, Pin+=NbytesIn, Pout+=NbytesOut) { @@ -296,10 +396,4 @@ } } - if (!strcmp (outtype, "byte") && !strcmp (intype, "byte")) { - int NbytesIn = 1; - for (i = 0; i < Nval*Nrow; i++, Pin+=NbytesIn, Pout+=NbytesOut) { - *(char *)Pout = (*(char *)Pin - Bzero) / Bscale; - } - } if (!strcmp (outtype, "short") && !strcmp (intype, "byte")) { int NbytesIn = 1; @@ -307,5 +401,5 @@ *(short *)Pout = (*(char *)Pin - Bzero) / Bscale; # ifdef BYTE_SWAP - SWAP_BYTE; + if (!nativeOrder) { SWAP_BYTE; } # endif } @@ -316,5 +410,5 @@ *(int *)Pout = (*(char *)Pin - Bzero) / Bscale; # ifdef BYTE_SWAP - SWAP_WORD; + if (!nativeOrder) { SWAP_WORD; } # endif } @@ -325,5 +419,5 @@ *(int64_t *)Pout = (*(char *)Pin - Bzero) / Bscale; # ifdef BYTE_SWAP - SWAP_DBLE; + if (!nativeOrder) { SWAP_DBLE; } # endif } @@ -334,5 +428,5 @@ *(float *)Pout = (*(char *)Pin - Bzero) / Bscale; # ifdef BYTE_SWAP - SWAP_WORD; + if (!nativeOrder) { SWAP_WORD; } # endif } @@ -343,5 +437,5 @@ *(double *)Pout = (*(char *)Pin - Bzero) / Bscale; # ifdef BYTE_SWAP - SWAP_DBLE; + if (!nativeOrder) { SWAP_DBLE; } # endif } @@ -366,5 +460,5 @@ *(short *)Pout = (*(short *)Pin - Bzero) / Bscale; # ifdef BYTE_SWAP - SWAP_BYTE; + if (!nativeOrder) { SWAP_BYTE; } # endif } @@ -375,5 +469,5 @@ *(int *)Pout = (*(short *)Pin - Bzero) / Bscale; # ifdef BYTE_SWAP - SWAP_WORD; + if (!nativeOrder) { SWAP_WORD; } # endif } @@ -384,5 +478,5 @@ *(int64_t *)Pout = (*(short *)Pin - Bzero) / Bscale; # ifdef BYTE_SWAP - SWAP_DBLE; + if (!nativeOrder) { SWAP_DBLE; } # endif } @@ -393,5 +487,5 @@ *(float *)Pout = (*(short *)Pin - Bzero) / Bscale; # ifdef BYTE_SWAP - SWAP_WORD; + if (!nativeOrder) { SWAP_WORD; } # endif } @@ -402,5 +496,5 @@ *(double *)Pout = (*(short *)Pin - Bzero) / Bscale; # ifdef BYTE_SWAP - SWAP_DBLE; + if (!nativeOrder) { SWAP_DBLE; } # endif } @@ -425,5 +519,5 @@ *(short *)Pout = (*(int *)Pin - Bzero) / Bscale; # ifdef BYTE_SWAP - SWAP_BYTE; + if (!nativeOrder) { SWAP_BYTE; } # endif } @@ -434,5 +528,5 @@ *(int *)Pout = (*(int *)Pin - Bzero) / Bscale; # ifdef BYTE_SWAP - SWAP_WORD; + if (!nativeOrder) { SWAP_WORD; } # endif } @@ -443,5 +537,5 @@ *(int64_t *)Pout = (*(int *)Pin - Bzero) / Bscale; # ifdef BYTE_SWAP - SWAP_DBLE; + if (!nativeOrder) { SWAP_DBLE; } # endif } @@ -452,5 +546,5 @@ *(float *)Pout = (*(int *)Pin - Bzero) / Bscale; # ifdef BYTE_SWAP - SWAP_WORD; + if (!nativeOrder) { SWAP_WORD; } # endif } @@ -461,5 +555,5 @@ *(double *)Pout = (*(int *)Pin - Bzero) / Bscale; # ifdef BYTE_SWAP - SWAP_DBLE; + if (!nativeOrder) { SWAP_DBLE; } # endif } @@ -484,5 +578,5 @@ *(short *)Pout = (*(int64_t *)Pin - Bzero) / Bscale; # ifdef BYTE_SWAP - SWAP_BYTE; + if (!nativeOrder) { SWAP_BYTE; } # endif } @@ -493,5 +587,5 @@ *(int *)Pout = (*(int64_t *)Pin - Bzero) / Bscale; # ifdef BYTE_SWAP - SWAP_WORD; + if (!nativeOrder) { SWAP_WORD; } # endif } @@ -502,5 +596,5 @@ *(int64_t *)Pout = (*(int64_t *)Pin - Bzero) / Bscale; # ifdef BYTE_SWAP - SWAP_DBLE; + if (!nativeOrder) { SWAP_DBLE; } # endif } @@ -511,5 +605,5 @@ *(float *)Pout = (*(int64_t *)Pin - Bzero) / Bscale; # ifdef BYTE_SWAP - SWAP_WORD; + if (!nativeOrder) { SWAP_WORD; } # endif } @@ -520,5 +614,5 @@ *(double *)Pout = (*(int64_t *)Pin - Bzero) / Bscale; # ifdef BYTE_SWAP - SWAP_DBLE; + if (!nativeOrder) { SWAP_DBLE; } # endif } @@ -543,5 +637,5 @@ *(short *)Pout = (*(float *)Pin - Bzero) / Bscale; # ifdef BYTE_SWAP - SWAP_BYTE; + if (!nativeOrder) { SWAP_BYTE; } # endif } @@ -552,5 +646,5 @@ *(int *)Pout = (*(float *)Pin - Bzero) / Bscale; # ifdef BYTE_SWAP - SWAP_WORD; + if (!nativeOrder) { SWAP_WORD; } # endif } @@ -561,5 +655,5 @@ *(int64_t *)Pout = (*(float *)Pin - Bzero) / Bscale; # ifdef BYTE_SWAP - SWAP_DBLE; + if (!nativeOrder) { SWAP_DBLE; } # endif } @@ -570,5 +664,5 @@ *(float *)Pout = (*(float *)Pin - Bzero) / Bscale; # ifdef BYTE_SWAP - SWAP_WORD; + if (!nativeOrder) { SWAP_WORD; } # endif } @@ -579,5 +673,5 @@ *(double *)Pout = (*(float *)Pin - Bzero) / Bscale; # ifdef BYTE_SWAP - SWAP_DBLE; + if (!nativeOrder) { SWAP_DBLE; } # endif } @@ -602,5 +696,5 @@ *(short *)Pout = (*(double *)Pin - Bzero) / Bscale; # ifdef BYTE_SWAP - SWAP_BYTE; + if (!nativeOrder) { SWAP_BYTE; } # endif } @@ -611,5 +705,5 @@ *(int *)Pout = (*(double *)Pin - Bzero) / Bscale; # ifdef BYTE_SWAP - SWAP_WORD; + if (!nativeOrder) { SWAP_WORD; } # endif } @@ -620,5 +714,5 @@ *(int64_t *)Pout = (*(double *)Pin - Bzero) / Bscale; # ifdef BYTE_SWAP - SWAP_DBLE; + if (!nativeOrder) { SWAP_DBLE; } # endif } @@ -629,5 +723,5 @@ *(float *)Pout = (*(double *)Pin - Bzero) / Bscale; # ifdef BYTE_SWAP - SWAP_WORD; + if (!nativeOrder) { SWAP_WORD; } # endif } @@ -638,8 +732,9 @@ *(double *)Pout = (*(double *)Pin - Bzero) / Bscale; # ifdef BYTE_SWAP - SWAP_DBLE; -# endif - } - } + if (!nativeOrder) { SWAP_DBLE; } +# endif + } + } +# endif /* check array space */ @@ -683,5 +778,5 @@ int gfits_set_table_column (Header *header, FTable *table, char *label, void *data, off_t Nrow) { - off_t Nx, Ny, nbytes; + off_t Nx, Ny; int i, N, Nfields, Nval, Nbytes, Nstart, Nv, Nb; char tlabel[80], field[80], format[80], cformat[80], type[16], tmp[16]; @@ -713,11 +808,5 @@ if (Ny == 0) { Ny = Nrow; - header[0].Naxis[1] = Ny; - gfits_modify (header, "NAXIS2", OFF_T_FMT, 1, Ny); - - nbytes = gfits_data_size (header); - REALLOCATE (table[0].buffer, char, MAX (nbytes, 1)); - bzero (table[0].buffer, nbytes); - table[0].datasize = nbytes; + gfits_set_table_rows (header, table, Ny); } if (Ny != Nrow) return (FALSE); Index: trunk/Ohana/src/libfits/table/F_table_format.c =================================================================== --- trunk/Ohana/src/libfits/table/F_table_format.c (revision 38255) +++ trunk/Ohana/src/libfits/table/F_table_format.c (revision 38441) @@ -8,13 +8,12 @@ int gfits_bintable_format (char *format, char *type, int *Nval, int *Nbytes) { - char *Lchar, *Fchar; - int Nv; - - if (format == (char *) NULL) return (FALSE); - if (format[0] == 0) return (FALSE); - Lchar = &format[strlen(format) - 1]; - - Nv = strtol (format, &Fchar, 10); - if (Fchar != Lchar) return (FALSE); + char *Fchar; + + if (format == (char *) NULL) return (FALSE); // format must be defined + if (format[0] == 0) return (FALSE); // format must not be empty string + // char *Lchar = &format[strlen(format) - 1]; // last-char pointer + + int Nv = strtol (format, &Fchar, 10); // Fchar points at end of leading number + // if (Fchar != Lchar) return (FALSE); -- this test is invalid for eg 1PB(nn) if ((Nv == 0) && (Fchar == format)) Nv = 1; @@ -51,7 +50,4 @@ { *Nbytes = 8; strcpy (type, "double"); *Nval = Nv; } break; - case 'P': - { *Nbytes = 8; strcpy (type, "var"); *Nval = 2*Nv; } - break; case 'C': { *Nbytes = 8; strcpy (type, "float"); *Nval = 2*Nv; } @@ -59,4 +55,10 @@ case 'M': { *Nbytes = 16; strcpy (type, "double"); *Nval = 2*Nv; } + break; + case 'P': + { *Nbytes = 4; strcpy (type, "var"); *Nval = 2*Nv; } + break; + case 'Q': + { *Nbytes = 8; strcpy (type, "var64"); *Nval = 2*Nv; } break; default: @@ -80,9 +82,9 @@ C - complex float M - complex double - P - var length array descrpt (64 bit) + P - var length array descrpt (32 bit pointer) + Q - var length array descrpt (64 bit pointer) all can be preceeded by integer N which specified number of entries in array - */ Index: trunk/Ohana/src/libfits/table/F_table_row.c =================================================================== --- trunk/Ohana/src/libfits/table/F_table_row.c (revision 38255) +++ trunk/Ohana/src/libfits/table/F_table_row.c (revision 38441) @@ -28,4 +28,5 @@ REALLOCATE (table[0].buffer, char, MAX (nbytes, 1)); table[0].datasize = nbytes; + table[0].heap_start = gfits_heap_start (header); memcpy (&table[0].buffer[Nstart], data, Nbytes*Nrow); @@ -67,4 +68,5 @@ table[0].datasize = gfits_data_size (table[0].header); + // table[0].heap_start = gfits_heap_start (header); table[0].pad = table[0].datasize - Nx*Ny; @@ -102,4 +104,5 @@ REALLOCATE (table[0].buffer, char, MAX (nbytes, 1)); table[0].datasize = nbytes; + table[0].heap_start = gfits_heap_start (header); return (TRUE); } Index: trunk/Ohana/src/libfits/table/F_table_varlength.c =================================================================== --- trunk/Ohana/src/libfits/table/F_table_varlength.c (revision 38255) +++ trunk/Ohana/src/libfits/table/F_table_varlength.c (revision 38441) @@ -2,4 +2,5 @@ # include +// fills the column definition structure for a specific column based on a given table. int gfits_varlength_column_define (FTable *ftable, VarLengthColumn *def, int column) { @@ -14,4 +15,6 @@ snprintf (field, 80, "TFORM%d", column); if (!gfits_scan (ftable->header, field, "%s", 1, format)) return (FALSE); + + def->column = column; // find and remove the max field length element @@ -48,4 +51,5 @@ if (*p1 == 'I') { def->nbytes = 2; def->format = *p1; } if (*p1 == 'J') { def->nbytes = 4; def->format = *p1; } + if (*p1 == 'K') { def->nbytes = 8; def->format = *p1; } if (*p1 == 'E') { def->nbytes = 4; def->format = *p1; } if (*p1 == 'D') { def->nbytes = 8; def->format = *p1; } @@ -54,18 +58,58 @@ if (!def->nbytes) { return (FALSE); } - /* scan columns to find column offset */ - def->Nstart = 0; + /* scan columns to find column offset for metadata column*/ + def->offset = 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->Nstart += Nv*Nb; + def->offset += Nv*Nb; } // heap_start must be long long so file may be very large - if (!gfits_scan (ftable->header, "THEAP", OFF_T_FMT, 1, &def->heap_start)) { - def->heap_start = ftable->header->Naxis[0]*ftable->header->Naxis[1]; - } - + // confirm that ftable->heap_start is correctly set? + + return TRUE; +} + +# define SWAP_WORD \ + tmp = pchar[0]; pchar[0] = pchar[3]; pchar[3] = tmp; \ + tmp = pchar[1]; pchar[1] = pchar[2]; pchar[2] = tmp; + +# define SWAP_DBLE \ + tmp = pchar[0]; pchar[0] = pchar[7]; pchar[7] = tmp; \ + tmp = pchar[1]; pchar[1] = pchar[6]; pchar[6] = tmp; \ + tmp = pchar[2]; pchar[2] = pchar[5]; pchar[5] = tmp; \ + tmp = pchar[3]; pchar[3] = pchar[4]; pchar[4] = tmp; + +int gfits_byteswap_varlength_column (FTable *ftable, int column) { + +# ifdef BYTE_SWAP + VarLengthColumn zdef; + if (!gfits_varlength_column_define (ftable, &zdef, column)) return FALSE; + + char *bufstart = &ftable->buffer[zdef.offset]; + + int i; + char *pchar, tmp; + + off_t Nx = ftable->header->Naxis[0]; + + for (i = 0; i < ftable->header->Naxis[1]; i++) { + if (zdef.mode == 'P') { + pchar = &bufstart[i*Nx]; + SWAP_WORD; + pchar = &bufstart[i*Nx + 4]; + SWAP_WORD; + } + if (zdef.mode == 'Q') { + char *pchar, tmp; + pchar = &bufstart[i*Nx]; + SWAP_DBLE; + pchar = &bufstart[i*Nx + 8]; + SWAP_DBLE; + } + } +# endif return TRUE; } @@ -75,41 +119,106 @@ void *gfits_varlength_column_pointer (FTable *ftable, VarLengthColumn *column, off_t row, off_t *length) { - void *result; - off_t offset, Nx; - if ((column->mode != 'P') && (column->mode != 'Q')) abort(); // find the values for the specified row // the values in the main table for this row and varlength column: - 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 + off_t Nx = ftable->header->Naxis[0]; + off_t offset = 0; if (column->mode == 'P') { int *ptr; - ptr = (int *) &ftable->buffer[row*Nx + column->Nstart]; + ptr = (int *) &ftable->buffer[row*Nx + column->offset]; *length = ptr[0]; offset = ptr[1]; - } else if (column->mode == 'Q') { + } + if (column->mode == 'Q') { off_t *ptr; - ptr = (off_t *) &ftable->buffer[row*Nx + column->Nstart]; + ptr = (off_t *) &ftable->buffer[row*Nx + column->offset]; *length = ptr[0]; offset = ptr[1]; - } else { - // this was actually trapped above - abort(); - } - - result = (void *) (ftable->buffer + column->heap_start + offset); + } + + // fprintf (stderr, "length: %d, offset: %d\n", (int) *length, (int) offset); + + void *result = (void *) (ftable->buffer + ftable->heap_start + offset); return result; } + +// to generate a compressed image or a compressed table, we need to be able to add data +// for varlength columns to the heap. To start, assume we have defined the cartesian +// portion of the table correctly, and are only extending the heap. + +int gfits_varlength_column_add_data (FTable *table, char *data, off_t Ndata, int row, VarLengthColumn *column) { + + // find the current starting point for new data (end of current buffer main data + current heap + off_t heap_offset = gfits_data_min_size (table->header); + + // extend the buffer size to add the new data + off_t Nbytes = gfits_data_pad_size (heap_offset + Ndata); + REALLOCATE (table->buffer, char, Nbytes); + + char *heap_ptr = table->buffer + heap_offset; + + // add the new data to the table buffer + memcpy (heap_ptr, data, Ndata); + + // *** now we add the data description to the table *** + + if ((column->mode != 'P') && (column->mode != 'Q')) abort(); + + // find the cell for the specified row & column + // the values in the main table for this row and varlength column: + int Nx = table->header->Naxis[0]; + + // set the cell values + if (column->mode == 'P') { + int *ptr; + ptr = (int *) &table->buffer[row*Nx + column->offset]; + ptr[0] = Ndata; + ptr[1] = heap_offset - table->heap_start; + } + if (column->mode == 'Q') { + off_t *ptr; + ptr = (off_t *) &table->buffer[row*Nx + column->offset]; + ptr[0] = Ndata; + ptr[1] = heap_offset - table->heap_start; + } + + // update maxlen + // fprintf (stderr, "row: %d, Ndata: %d, maxlen: %d\n", (int) row, (int) Ndata, (int) column->maxlen); + column->maxlen = MAX (column->maxlen, Ndata); + + off_t pcount; + gfits_scan (table->header, "PCOUNT", OFF_T_FMT, 1, &pcount); + + pcount += Ndata; + gfits_modify (table->header, "PCOUNT", OFF_T_FMT, 1, pcount); + table->header->pcount = pcount; + table->datasize = gfits_data_size (table->header); + myAssert (table->datasize == Nbytes, "inconsistent header and data"); + + return TRUE; +} + +// update the header TFORM value for this entry +int gfits_varlength_column_finish (FTable *ftable, VarLengthColumn *def) { + + char field[81]; + char format[81]; + + // grab the value of TFORMn for this column + snprintf (field, 80, "TFORM%d", def->column); + if (!gfits_scan (ftable->header, field, "%s", 1, format)) return FALSE; + + // find and update the max field length element + char *p1 = strchr (format, '('); + if (!p1) return (FALSE); // not a valid varlength column -- missing (e_max) + + p1 ++; + *p1 = 0; + + snprintf (p1, 80 - (p1 - format), "%d)", def->maxlen); + if (!gfits_modify (ftable->header, field, "%s", 1, format)) return FALSE; + + return TRUE; +} Index: trunk/Ohana/src/libfits/table/F_uncompress_T.c =================================================================== --- trunk/Ohana/src/libfits/table/F_uncompress_T.c (revision 38441) +++ trunk/Ohana/src/libfits/table/F_uncompress_T.c (revision 38441) @@ -0,0 +1,424 @@ +# include +# include +# include + +# define ESCAPE { fprintf (stderr, "error in %s @ line %d\n", __func__, __LINE__); goto escape; } + +int gfits_distribute_table_data (FTable *table, TableField *field, char *raw, int row_start, int Nrows); +int gfits_distribute_table_gzp2 (FTable *table, TableField *field, char *raw, int row_start, int Nrows); + +# define VERBOSE_DUMP 0 + +int gfits_dump_raw_table (FTable *table, char *message) { +# if (VERBOSE_DUMP) + int i; + fprintf (stderr, "%s data:\n", message); + for (i = 0; i < table->header->Naxis[0]*table->header->Naxis[1]; i++) { + // for (i = 0; i < 16; i++) { + fprintf (stderr, "%02hhx", table->buffer[i]); + if (i % 2) fprintf (stderr, " "); + if (i % 32 == 31) fprintf (stderr, "\n"); + } + fprintf (stderr, "\n"); +# endif + return TRUE; +} + +int gfits_dump_cmp_table (FTable *table, char *message) { +# if (VERBOSE_DUMP) + int i; + fprintf (stderr, "%s pntr:\n", message); + for (i = 0; i < table->header->Naxis[0]*table->header->Naxis[1]; i++) { + fprintf (stderr, "%02hhx", table->buffer[i]); + if (i % 2) fprintf (stderr, " "); + if (i % 32 == 31) fprintf (stderr, "\n"); + } + fprintf (stderr, "\n"); + + fprintf (stderr, "%s data:\n", message); + for (i = 0; i < table->header->pcount; i++) { + fprintf (stderr, "%02hhx", table->buffer[table->heap_start + i]); + if (i % 2) fprintf (stderr, " "); + if (i % 32 == 31) fprintf (stderr, "\n"); + } + fprintf (stderr, "\n"); +# endif + return TRUE; +} + +static float timeSum1 = 0.0; +static float timeSum2 = 0.0; +static float timeSum2a = 0.0; +static float timeSum2b = 0.0; +static float timeSum2c = 0.0; +static float timeSum2d = 0.0; +static float timeSum2e = 0.0; +static float timeSum2f = 0.0; +static float timeSum3 = 0.0; + +void gfits_uncompress_timing () { + fprintf (stderr, "unc times: %f %f %f\n", timeSum1, timeSum2, timeSum3); + fprintf (stderr, "unc times: %f %f %f %f %f %f\n", timeSum2a, timeSum2b, timeSum2c, timeSum2d, timeSum2e, timeSum2f); + + timeSum1 = 0.0; + timeSum2 = 0.0; + timeSum2a = 0.0; + timeSum2b = 0.0; + timeSum2c = 0.0; + timeSum2d = 0.0; + timeSum2e = 0.0; + timeSum2f = 0.0; + timeSum3 = 0.0; +} + +int gfits_uncompress_table (FTable *srctable, FTable *tgttable) { + + int i; + + char keyword[81]; + + char *raw = NULL; + char *zdata = NULL; + + Header *srcheader = srctable->header; + Header *tgtheader = tgttable->header; + + // XXX EAM: + struct timeval startTimer, stopTimer; + float dtime; + gettimeofday (&startTimer, (void *) NULL); + + int ztilelen; + if (!gfits_scan (srcheader, "ZTILELEN", "%d", 1, &ztilelen)) ESCAPE; + // XXX if ZTILELEN is missing, it should have the value of ZNAXIS2 + + // EXTNAME is not mandatory, keep it if it exists: + char extname[80]; + if (!gfits_scan (srcheader, "EXTNAME", "%s", 1, extname)) { + strcpy (extname, "UNCOMPRESSED_TABLE"); + } + + // creates an empty (Naxes = 2, Naxis[i] = 0) table header with XTENSION = BINTABLE, EXTNAME = COMPRESSED_TABLE + if (!gfits_create_table_header (tgtheader, "BINTABLE", extname)) ESCAPE; + + int Nfields; + if (!gfits_scan (srcheader, "TFIELDS", "%d", 1, &Nfields)) ESCAPE; + + // below we generate a fake table with tmpbuffers. allocated it to the max needed size + int max_width = 0; + int offset = 0; + + TableField *fields = NULL; + ALLOCATE (fields, TableField, Nfields); + for (i = 0; i < Nfields; i++) { + snprintf (keyword, 80, "TTYPE%d", i+1); + if (!gfits_scan (srcheader, keyword, "%s", 1, fields[i].ttype)) ESCAPE; + if (!gfits_scan_alt (srcheader, keyword, "%C", 1, fields[i].ttype_cmt)) ESCAPE; + + // TUNIT is not mandatory + snprintf (keyword, 80, "TUNIT%d", i+1); + if (!gfits_scan (srcheader, keyword, "%s", 1, fields[i].tunit)) { + fields[i].tunit[0] = 0; + } else { + if (!gfits_scan_alt (srcheader, keyword, "%C", 1, fields[i].tunit_cmt)) ESCAPE; + } + + snprintf (keyword, 80, "ZFORM%d", i+1); + if (!gfits_scan (srcheader, keyword, "%s", 1, fields[i].tformat)) ESCAPE; + if (!gfits_scan_alt (srcheader, keyword, "%C", 1, fields[i].tformat_cmt)) ESCAPE; + + snprintf (keyword, 80, "ZCTYP%d", i+1); + if (!gfits_scan (srcheader, keyword, "%s", 1, fields[i].zctype)) { + strcpy (fields[i].zctype, "GZIP_2"); // XXX NOTE: not yet supported + } + + if (!gfits_bintable_format (fields[i].tformat, fields[i].datatype, &fields[i].Nvalues, &fields[i].pixsize)) ESCAPE; + fields[i].rowsize = fields[i].Nvalues*fields[i].pixsize; + max_width = MAX(max_width, fields[i].rowsize); + + fields[i].offset = offset; + offset += fields[i].rowsize; + + if (!gfits_varlength_column_define (srctable, &fields[i].zdef, i + 1)) ESCAPE; + + if (!gfits_define_bintable_column (tgtheader, fields[i].tformat, fields[i].ttype, fields[i].ttype_cmt, fields[i].tunit, 1.0, 0.0)) ESCAPE; + } + + // allocates the default data array (tgttable->buffer) + if (!gfits_create_table (tgtheader, tgttable)) ESCAPE; + + // copy original header to output header (XXX this needs to be finished) + if (!gfits_copy_keywords_compress (srcheader, tgtheader)) ESCAPE; + + off_t Nx, Ny; + if (!gfits_scan (srcheader, "ZNAXIS1", OFF_T_FMT, 1, &Nx)) ESCAPE; + if (!gfits_scan (srcheader, "ZNAXIS2", OFF_T_FMT, 1, &Ny)) ESCAPE; + myAssert (Nx == tgtheader->Naxis[0], "table definition error?"); + + // XXX TIMER 1 + gettimeofday (&stopTimer, (void *) NULL); + dtime = DTIME (stopTimer, startTimer); + timeSum1 += dtime; + gettimeofday (&startTimer, (void *) NULL); + + // this allocates the pointer data array for the full set of tiles (NOT the heap) + char *tmpbuffer = NULL; + ALLOCATE_ZERO (tmpbuffer, char, max_width*Ny); + for (i = 0; i < Nfields; i++) { + if (!gfits_set_bintable_column (tgtheader, tgttable, fields[i].ttype, tmpbuffer, Ny)) ESCAPE; + } + free (tmpbuffer); + + int Ntile = srcheader->Naxis[1]; + int ztilelast; + if (Ntile == 0) { + ztilelast = 0; + } else { + ztilelast = (tgtheader->Naxis[1] % ztilelen) ? tgtheader->Naxis[1] % ztilelen : ztilelen; + } + + // allocate the intermediate storage buffers + int Nraw_alloc = max_width*ztilelen + 1000; + ALLOCATE (raw, char, Nraw_alloc); + ALLOCATE (zdata, char, max_width*ztilelen); + + // uncompress the data : extract a tile's compressed data, uncompress the tile, insert into full table + // each tile -> 1 row of the output table + + gfits_dump_cmp_table (srctable, "unc"); + + // XXX TIMER 2 + gettimeofday (&stopTimer, (void *) NULL); + dtime = DTIME (stopTimer, startTimer); + timeSum2 += dtime; + gettimeofday (&startTimer, (void *) NULL); + + off_t row; + for (row = 0; row < srcheader->Naxis[1]; row++) { + + for (i = 0; i < Nfields; i++) { + + // gfits_uncompress_data can take values specific to the compression mode + // optname, optvalue = NULL, Noptions = 0 + + // XXX note that this function currently byteswaps. I should redo the code to make + // byteswap of the varlength data segment a separate function + + gettimeofday (&startTimer, (void *) NULL); + + off_t Nzdata; + char *zdata = gfits_varlength_column_pointer (srctable, &fields[i].zdef, row, &Nzdata); + if (!zdata) ESCAPE; + + if (VERBOSE_DUMP) { + int k; + fprintf (stderr, "u4: "); + for (k = 0; k < Nzdata; k++) { + fprintf (stderr, "%02hhx", zdata[k]); + if (k % 2) fprintf (stderr, " "); + // if (k % 32 == 31) fprintf (stderr, "\n"); + } + fprintf (stderr, "\n"); + } + + // XXX TIMER 2a + gettimeofday (&stopTimer, (void *) NULL); + dtime = DTIME (stopTimer, startTimer); + timeSum2a += dtime; + gettimeofday (&startTimer, (void *) NULL); + + if (strcasecmp(fields[i].zctype, "NONE_2") && + strcasecmp(fields[i].zctype, "GZIP_1") && + strcasecmp(fields[i].zctype, "GZIP_2") && + strcasecmp(fields[i].zctype, "RICE_1") && + strcasecmp(fields[i].zctype, "RICE_ONE")) { + // Nzdata is number of bytes + if (!gfits_byteswap_zdata (zdata, Nzdata, fields[i].pixsize)) ESCAPE; + } + + // XXX TIMER 2b + gettimeofday (&stopTimer, (void *) NULL); + dtime = DTIME (stopTimer, startTimer); + timeSum2b += dtime; + gettimeofday (&startTimer, (void *) NULL); + + if (VERBOSE_DUMP) { + int k; + fprintf (stderr, "u3: "); + for (k = 0; k < Nzdata; k++) { + fprintf (stderr, "%02hhx", zdata[k]); + if (k % 2) fprintf (stderr, " "); + // if (k % 32 == 31) fprintf (stderr, "\n"); + } + fprintf (stderr, "\n"); + } + + // XXX TIMER 2c + gettimeofday (&stopTimer, (void *) NULL); + dtime = DTIME (stopTimer, startTimer); + timeSum2c += dtime; + gettimeofday (&startTimer, (void *) NULL); + + int Nrows = (row == Ntile - 1) ? ztilelast : ztilelen; + int Nraw = Nrows*fields[i].Nvalues*fields[i].pixsize; // expected number of bytes + if (!gfits_uncompress_data (zdata, Nzdata, fields[i].zctype, NULL, NULL, 0, raw, &Nraw, Nraw_alloc, fields[i].pixsize)) ESCAPE; + + if (VERBOSE_DUMP) { + int k; + fprintf (stderr, "u2: "); + for (k = 0; k < Nraw*fields[i].pixsize; k++) { + fprintf (stderr, "%02hhx", raw[k]); + if (k % 2) fprintf (stderr, " "); + // if (k % 32 == 31) fprintf (stderr, "\n"); + } + fprintf (stderr, "\n"); + } + + // XXX TIMER 2d + gettimeofday (&stopTimer, (void *) NULL); + dtime = DTIME (stopTimer, startTimer); + timeSum2d += dtime; + gettimeofday (&startTimer, (void *) NULL); + + if (!strcasecmp(fields[i].zctype, "GZIP_1")) { + myAssert ((fields[i].zdef.format != 'C') && (fields[i].zdef.format != 'M'), "swap is probably wrong for C or M columns"); + if (!gfits_byteswap_zdata (raw, Nraw * fields[i].pixsize, fields[i].pixsize)) ESCAPE; + } + + // XXX TIMER 2e + gettimeofday (&stopTimer, (void *) NULL); + dtime = DTIME (stopTimer, startTimer); + timeSum2e += dtime; + gettimeofday (&startTimer, (void *) NULL); + + if (VERBOSE_DUMP) { + int k; + fprintf (stderr, "u1: "); + for (k = 0; k < Nraw*fields[i].pixsize; k++) { + fprintf (stderr, "%02hhx", raw[k]); + if (k % 2) fprintf (stderr, " "); + // if (k % 32 == 31) fprintf (stderr, "\n"); + } + fprintf (stderr, "\n"); + } + + // int valid = (row == 0) && (i == 0); + + off_t row_start = row*ztilelen; + // ^- first row for this tile & field + // NOTE: assumes each tile has the same length, ex the last (true for now) + + // copy the raw pixels from their native matrix locations to the temporary output buffer + if (!strcasecmp(fields[i].zctype, "GZIP_2") || !strcasecmp(fields[i].zctype, "NONE_2")) { + if (!gfits_distribute_table_gzp2 (tgttable, &fields[i], raw, row_start, Nrows)) ESCAPE; + } else { + if (!gfits_distribute_table_data (tgttable, &fields[i], raw, row_start, Nrows)) ESCAPE; + } + + // XXX TIMER 2f + gettimeofday (&stopTimer, (void *) NULL); + dtime = DTIME (stopTimer, startTimer); + timeSum2f += dtime; + gettimeofday (&startTimer, (void *) NULL); + } + } + + // XXX TIMER 3 + gettimeofday (&stopTimer, (void *) NULL); + dtime = DTIME (stopTimer, startTimer); + timeSum3 += dtime; + gettimeofday (&startTimer, (void *) NULL); + + gfits_dump_raw_table (tgttable, "unc"); + + free (raw); + free (zdata); + return (TRUE); + +escape: + FREE (raw); + FREE (zdata); + return FALSE; +} + +# define VERBOSE 0 + +// 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_distribute_table_data (FTable *table, TableField *field, char *raw, int row_start, int Nrows) { + + off_t i; + + // we are copying NN rows into the column which starts at XX and has MM bytes per row + + off_t Nx = table->header->Naxis[0]; + + char *tblbuffer = &table->buffer[Nx*row_start + field->offset]; + + int rowsize = field->rowsize; + + if (VERBOSE) fprintf (stderr, "distribute: "); + for (i = 0; i < Nrows; i++, tblbuffer += Nx) { + memcpy (tblbuffer, &raw[i*rowsize], rowsize); +# if (VERBOSE) + int j; for (j = 0; j < field->rowsize; j++) fprintf (stderr, "0x%02hhx ", table->buffer[Nx*row + field->offset + j]); +# endif + } + if (VERBOSE) fprintf (stderr, "\n"); + + return (TRUE); +} + +int gfits_distribute_table_gzp2 (FTable *table, TableField *field, char *raw, int row_start, int Nrows) { + + off_t i, j, k; + + // we are copying NN rows into the column which starts at XX and has MM bytes per row + + off_t Nx = table->header->Naxis[0]; + int Nvalues = field->Nvalues; + int pixsize = field->pixsize; + int offset = field->offset; + + char *rawptr = raw; + + for (k = 0; k < field->pixsize; k++) { +# ifdef BYTE_SWAP + char *tblptr_start = &table->buffer[Nx*row_start + offset + (pixsize - k - 1)]; +# else + char *tblptr_start = &table->buffer[Nx*row_start + offset + k]; +# endif + for (i = 0; i < Nrows; i++, tblptr_start += Nx) { + char *tblptr = tblptr_start; + for (j = 0; j < Nvalues; j++, tblptr += pixsize, rawptr++) { + *tblptr = *rawptr; + } + } + } + return (TRUE); +} + +int gfits_distribute_table_gzp2_alt (FTable *table, TableField *field, char *raw, int row_start, int Nrows) { + + off_t i, j, k; + + // we are copying NN rows into the column which starts at XX and has MM bytes per row + + off_t Nx = table->header->Naxis[0]; + + for (k = 0; k < field->pixsize; k++) { + for (i = 0; i < Nrows; i++) { + int row = row_start + i; + char *rawptr = &raw[i*field->Nvalues + k*Nrows*field->Nvalues]; +# ifdef BYTE_SWAP + char *tblptr = &table->buffer[Nx*row + field->offset + (field->pixsize - k - 1)]; +# else + char *tblptr = &table->buffer[Nx*row + field->offset + k]; +# endif + for (j = 0; j < field->Nvalues; j++, tblptr += field->pixsize, rawptr++) { + *tblptr = *rawptr; + } + } + } + return (TRUE); +} Index: trunk/Ohana/src/libfits/table/F_write_T.c =================================================================== --- trunk/Ohana/src/libfits/table/F_write_T.c (revision 38255) +++ trunk/Ohana/src/libfits/table/F_write_T.c (revision 38441) @@ -14,4 +14,5 @@ status = fseeko (f, 0LL, SEEK_END); /* write table to end of file! */ + if (status) { perror ("fseeko: "); return (FALSE); } status = gfits_fwrite_table (f, table); @@ -26,6 +27,5 @@ Nbytes = fwrite (table[0].buffer, sizeof(char), table[0].datasize, f); - - if (Nbytes != table[0].datasize) return (FALSE); + if (Nbytes != table[0].datasize) { perror ("fwrite: "); return (FALSE); } return (TRUE); } @@ -45,4 +45,5 @@ /* file pointer is at beginning of desired table data */ start = ftello (f); + if (start < 0) { perror ("ftello: "); return FALSE; } for (i = 0; i < Nrow; i++) { @@ -50,5 +51,5 @@ fseeko (f, offset, SEEK_SET); Nbytes = fwrite (table[0].buffer[i], sizeof (char), Nx, f); - if (Nbytes != Nx) { return (FALSE); } + if (Nbytes != Nx) { perror ("fwrite: "); return (FALSE); } } @@ -58,7 +59,9 @@ offset = start + Nx*Ny; - fseeko (f, offset, SEEK_SET); + int status = fseeko (f, offset, SEEK_SET); + if (status) { perror ("fseeko: "); return FALSE; } + Nbytes = fwrite (pad, sizeof (char), Npad, f); - if (Nbytes != Npad) { return (FALSE); } + if (Nbytes != Npad) { perror ("fwrite: "); return (FALSE); } free (pad); @@ -84,4 +87,5 @@ Nx = ftable[0].header[0].Naxis[0]; // final output table size on disk ftable[0].datasize = gfits_data_size (ftable[0].header); + ftable[0].heap_start = gfits_heap_start (ftable[0].header); Nskip = start * Nx; @@ -96,4 +100,5 @@ // cursor must be at start of the table (after table header) if (fseeko (f, Nskip, SEEK_CUR)) { + perror ("fseeko: "); fprintf (stderr, "can't seek table start\n"); return (FALSE); @@ -101,7 +106,5 @@ Nwrite = fwrite (ftable[0].buffer, sizeof (char), Nbytes, f); - if (Nwrite != Nbytes) { - return (FALSE); - } + if (Nwrite != Nbytes) { perror ("fwrite: "); return (FALSE); } if (Ntotal >= Ndisk) { @@ -112,5 +115,5 @@ free (pad); - if (Nbytes != Npad) return (FALSE); + if (Nbytes != Npad) { perror ("fwrite: "); return (FALSE); } } Index: trunk/Ohana/src/libfits/table/F_write_TH.c =================================================================== --- trunk/Ohana/src/libfits/table/F_write_TH.c (revision 38255) +++ trunk/Ohana/src/libfits/table/F_write_TH.c (revision 38441) @@ -14,4 +14,5 @@ status = fseeko (f, 0LL, SEEK_END); /* write header to end of file! */ + if (status) { perror ("fseeko: "); return (FALSE); } status = gfits_fwrite_Theader (f, header); @@ -26,6 +27,6 @@ Nbytes = fwrite (header[0].buffer, sizeof(char), header[0].datasize, f); + if (Nbytes != header[0].datasize) { perror ("fwrite: "); return (FALSE); } - if (Nbytes != header[0].datasize) return (FALSE); return (TRUE); } Index: trunk/Ohana/src/libfits/test/compress.sh =================================================================== --- trunk/Ohana/src/libfits/test/compress.sh (revision 38441) +++ trunk/Ohana/src/libfits/test/compress.sh (revision 38441) @@ -0,0 +1,446 @@ +input tap.dvo + +macro tests + test_image_all + test_table_all_rnd +end + +macro test_image_all + + tapPLAN 88 + + foreach mode NONE NONE_1 NONE_2 GZIP_1 GZIP_2 RICE_1 RICE_ONE + # foreach mode GZIP_1 + # foreach mode GZIP_2 + # foreach mode RICE_1 + foreach bitpix 8 16 32 -32 -64 + # foreach bitpix 32 + if (($mode == RICE_1) && ($bitpix < 0)) continue ; # RICE_1 is for int only + if (($mode == RICE_ONE) && ($bitpix < 0)) continue ; # RICE_1 is for int only + tapDIAG 2 "===== bitpix = $bitpix, mode = $mode =====" + test_image_single $bitpix $mode + end + end + + tapDONE +end + +macro test_image_single + if ($0 != 3) + echo "USAGE: test_image_single (bitpix) (mode)" + break + end + + local bitpix cmpmode + $bitpix = $1 + $cmpmode = $2 + + delete -q x z b c d e + + mcreate z 8 8 + set x = xramp(z) + wd x test.raw.fits -bitpix $bitpix -bzero 0 -bscale 1 + wd x test.cmp.fits -bitpix $bitpix -bzero 0 -bscale 1 -compress-mode $cmpmode + rd b test.cmp.fits -x 0 + set d = x - b + stat -q d + tapOK {abs($MEAN) < 0.001} "read our cmp $cmpmode $bitpix (MEAN = $MEAN)" + tapOK {$SIGMA < 0.001} "read our cmp $cmpmode $bitpix (SIGMA = $SIGMA)" + + if ($cmpmode == NONE) return + if ($cmpmode == NONE_1) return + if ($cmpmode == NONE_2) return + + 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 -q d + tapOK {abs($MEAN) < 0.001} "read fpack $cmpmode $bitpix (MEAN = $MEAN)" + tapOK {$SIGMA < 0.001} "read fpack $cmpmode $bitpix (SIGMA = $SIGMA)" + + if ($cmpmode == RICE_1_0) + echo "*** funpack fails on RICE_1 ****" + return + 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 -q d + tapOK {abs($MEAN) < 0.001} "funpack our cmp $cmpmode $bitpix (MEAN = $MEAN)" + tapOK {$SIGMA < 0.001} "funpack our cmp $cmpmode $bitpix (SIGMA = $SIGMA)" +end + +macro test_table_single_mode + if ($0 != 2) + echo "USAGE: test_table_single_mode (mode)" + break + end + + local cmpmode + $cmpmode = $1 + + delete -q x y ID + + $NPT = 10 + create x 0 $NPT + set y = x^2 + create ID 1 {$NPT + 1} -int + + $fields = ID x y + $format = IEE + + write -fits test test.raw.tbl $fields -format $format + write -fits test test.cmp.tbl $fields -format $format -compress-mode $cmpmode + + foreach f $fields + set $f\_raw = $f + delete $f + end + + echo "===== raw =====" + data test.raw.tbl + read -fits test $fields + + foreach f $fields + set dv = $f - $f\_raw + vstat dv + end + delete -q $fields + + echo "===== cmp =====" + data test.cmp.tbl + read -fits test $fields + + foreach f $fields + set dv = $f - $f\_raw + vstat dv + end +end + +macro test_table_all_single_column + + tapPLAN 88 + + local format cmpmode + + # foreach cmpmode NONE GZIP_1 GZIP_2 RICE_1 RICE_ONE + foreach cmpmode GZIP_2 + foreach format B I J K E D + if (($cmpmode == RICE_1) && ($format == E)) continue ; # RICE_1 is for int only + if (($cmpmode == RICE_ONE) && ($format == E)) continue ; # RICE_1 is for int only + if (($cmpmode == RICE_1) && ($format == D)) continue ; # RICE_1 is for int only + if (($cmpmode == RICE_ONE) && ($format == D)) continue ; # RICE_1 is for int only + if (($cmpmode == RICE_1) && ($format == K)) continue ; # RICE_1 is for int only + if (($cmpmode == RICE_ONE) && ($format == K)) continue ; # RICE_1 is for int only + test_table_single_column $cmpmode $format + end + end + tapDONE +end + +macro test_table_single_column + if ($0 != 3) + echo "USAGE: test_table_single_column (mode) (format)" + break + end + + local cmpmode format + $cmpmode = $1 + $format = $2 + + tapDIAG 2 "===== testing $cmpmode $format =====" + + delete -q x dv + + $NPT = 10000 + + if ($format == B) + create x -128 128 + else + create x 0 $NPT + end + + $fields = x + + write -fits test test.raw.tbl $fields -format $format + write -fits test test.cmp.tbl $fields -format $format -compress-mode $cmpmode + + foreach f $fields + set $f\_raw = $f + delete $f + end + + tapDIAG 2 "===== $cmpmode $format raw =====" + data test.raw.tbl + read -fits test $fields + + foreach f $fields + set dv = $f - $f\_raw + vstat -q dv + tapOK {abs($MEAN) < 0.001} "read our raw $cmpmode $format (MEAN = $MEAN)" + tapOK {$SIGMA < 0.001} "read our raw $cmpmode $format (SIGMA = $SIGMA)" + end + delete -q $fields + + tapDIAG 2 "===== $cmpmode $format cmp =====" + data test.cmp.tbl + read -fits test $fields + + foreach f $fields + set dv = $f - $f\_raw + vstat -q dv + tapOK {abs($MEAN) < 0.001} "read our cmp $cmpmode $format (MEAN = $MEAN)" + tapOK {$SIGMA < 0.001} "read our cmp $cmpmode $format (SIGMA = $SIGMA)" + end + + if ($cmpmode == NONE) return; # funpack will not recognize NONE + if ($cmpmode == RICE_ONE) return; # funpack will not recognize RICE_ONE for tables + + tapDIAG 2 "===== $cmpmode $format fpack =====" + + # try to run fpack on ours + if (1) + tapEXEC /bin/cp -f test.raw.tbl test.fpk.tbl + tapEXEC fpack -table -F test.fpk.tbl + + data test.fpk.tbl + read -fits test $fields + + foreach f $fields + set dv = $f - $f\_raw + vstat -q dv + tapOK {abs($MEAN) < 0.001} "read fpack $cmpmode $format (MEAN = $MEAN)" + tapOK {$SIGMA < 0.001} "read fpack $cmpmode $format (SIGMA = $SIGMA)" + end + delete -q $fields + end + + break -auto off + echo $cmpmode $format + exec ls test.cmp.tbl test.fpk.tbl | fields -x 0 PCOUNT ZCTYP1 + break -auto on + + tapDIAG 2 "===== $cmpmode $format funpack =====" + + # try to run fpack on ours + tapEXEC /bin/cp -f test.cmp.tbl test.fun.tbl + tapEXEC funpack -F test.fun.tbl + + data test.fun.tbl + read -fits test $fields + + foreach f $fields + set dv = $f - $f\_raw + vstat -q dv + tapOK {abs($MEAN) < 0.001} "read funpack $cmpmode $format (MEAN = $MEAN)" + tapOK {$SIGMA < 0.001} "read funpack $cmpmode $format (SIGMA = $SIGMA)" + end + delete -q $fields +end + +macro test_table_all_multi_column + + local format cmpmode + + foreach cmpmode NONE GZIP_1 + foreach format B I J K E D + test_table_multi_column $cmpmode $format + end + end +end + +macro test_table_multi_column + if ($0 != 3) + echo "USAGE: test_table_multi_column (mode) (format)" + break + end + + local cmpmode + $cmpmode = $1 + + delete -q x y z + + $NPT = 10000 + create x 0 $NPT + set y = x^2 + set z = 0.01*x + + local format + + $fields = x y z + $format = $2\$2\$2 + + $fields = x y + $format = $2\$2 + + write -fits test test.raw.tbl $fields -format $format + write -fits test test.cmp.tbl $fields -format $format -compress-mode $cmpmode + + foreach f $fields + set $f\_raw = $f + delete $f + end + + echo "===== $cmpmode $format raw =====" + data test.raw.tbl + read -fits test $fields + + foreach f $fields + set dv = $f - $f\_raw + vstat dv + end + delete -q $fields + + echo "===== $cmpmode $format cmp =====" + data test.cmp.tbl + read -fits test $fields + + foreach f $fields + set dv = $f - $f\_raw + vstat dv + end +end + +macro test_table_all_rnd + + local format cmpmode + + foreach cmpmode NONE NONE_1 NONE_2 GZIP_1 GZIP_2 RICE_1 RICE_ONE + test_table_rnd $cmpmode + end +end + +macro test_table_rnd + if ($0 != 2) + echo "USAGE: test_table_rnd (mode)" + break + end + + local cmpmode + $cmpmode = $1 + + tapDIAG 2 "===== testing $cmpmode =====" + + delete -q Achar Ashort Aint Along + delete -q Ffloat Fdouble + + # lrnd generates rnd int between 0 and 0x7fffff (2^31 - 1 inclusive) + + $NPT = 1000000 + # $NPT = 30 + create ni 0 $NPT -int + create nf 0 $NPT + + # set Achar = lrnd(ni) & 0x7f + set Achar = int(255.999*(drnd(ni) - 0.5)) + set Ashort = lrnd(ni) & 0x7fff + set Aint = lrnd(ni) + + # I cannot represent the full dynamic range of 64 bit with mana vector ints, but + # I can with the (double) floating point vector. however, I get the wrong bits + set Along_dbl = 9.223371e18 * (2.0*drnd(nf) - 1.0); + + cast Along = Along_dbl as int64 + + set Ffloat_dbl = 3.402822e+38 * (2.0*drnd(nf) - 1.0) + cast Ffloat = Ffloat_dbl as single + set Fdouble = 1.797692e+308 * (2.0*drnd(nf) - 1.0) + + local format + + $fields = Achar Ashort Aint Along Ffloat Fdouble + $format = BIJKED + + #$fields = Along + #$format = K + + write -fits test test.raw.tbl $fields -format $format + write -fits test test.cmp.tbl $fields -format $format -compress-mode $cmpmode + + foreach f $fields + set $f\_raw = $f + delete $f + end + + tapDIAG 2 "===== $cmpmode $format raw =====" + data test.raw.tbl + read -fits test $fields + + foreach f $fields + set dv = $f - $f\_raw + vstat -q dv + tapOK {abs($MEAN) < 0.001} "read our raw $f $cmpmode $format (MEAN = $MEAN)" + tapOK {$SIGMA < 0.001} "read our raw $f $cmpmode $format (SIGMA = $SIGMA)" + end + delete -q $fields + + tapDIAG 2 "===== $cmpmode $format cmp =====" + data test.cmp.tbl + read -fits test $fields + + foreach f $fields + set dv = $f - $f\_raw + vstat -q dv + tapOK {abs($MEAN) < 0.001} "read our cmp $f $cmpmode $format (MEAN = $MEAN)" + tapOK {$SIGMA < 0.001} "read our cmp $f $cmpmode $format (SIGMA = $SIGMA)" + end + delete -q $fields + + tapDIAG 2 "===== $cmpmode $format fpack =====" + + # try to run fpack on ours + if (1) + tapEXEC /bin/cp -f test.raw.tbl test.fpk.tbl + tapEXEC fpack -table -F test.fpk.tbl + + data test.fpk.tbl + read -fits test $fields + + foreach f $fields + set dv = $f - $f\_raw + vstat -q dv + tapOK {abs($MEAN) < 0.001} "read fpack $f $cmpmode $format (MEAN = $MEAN)" + tapOK {$SIGMA < 0.001} "read fpack $f $cmpmode $format (SIGMA = $SIGMA)" + end + delete -q $fields + end + + break -auto off + echo $cmpmode $format + exec ls test.cmp.tbl test.fpk.tbl | fields -x 0 PCOUNT ZCTYP1 ZCTYP2 ZCTYP3 ZCTYP4 ZCTYP5 ZCTYP6 + break -auto on + + if ($cmpmode == NONE) return; # funpack will not recognize NONE + if ($cmpmode == NONE_1) return; # funpack will not recognize NONE + if ($cmpmode == NONE_2) return; # funpack will not recognize NONE + if ($cmpmode == RICE_ONE) return; # funpack will not recognize RICE_ONE for tables + + tapDIAG 2 "===== $cmpmode $format funpack =====" + + # try to run fpack on ours + tapEXEC /bin/cp -f test.cmp.tbl test.fun.tbl + tapEXEC funpack -F test.fun.tbl + + data test.fun.tbl + read -fits test $fields + + foreach f $fields + set dv = $f - $f\_raw + vstat -q dv + tapOK {abs($MEAN) < 0.001} "read funpack $f $cmpmode $format (MEAN = $MEAN)" + tapOK {$SIGMA < 0.001} "read funpack $f $cmpmode $format (SIGMA = $SIGMA)" + end + delete -q $fields +end + Index: trunk/Ohana/src/libfits/test/imagecomp.c =================================================================== --- trunk/Ohana/src/libfits/test/imagecomp.c (revision 38441) +++ trunk/Ohana/src/libfits/test/imagecomp.c (revision 38441) @@ -0,0 +1,325 @@ +# include +# include +# include "tap_ohana.h" + +int test_compress (int bitpix, char *zcmptype); +int test_compress_fullrange (int bitpix, char *zcmptype); +int test_compress_fulltile (int bitpix, char *zcmptype); + +// char *cmptype[] = {"NONE", "NONE_2", "GZIP_1", "GZIP_2", "PLIO_1", "RICE_1", "RICE_ONE", "HCOMPRESS_1", NULL}; +char *cmptype[] = {"NONE", "NONE_1", "NONE_2", "GZIP_1", "GZIP_2", "RICE_1", "RICE_ONE", NULL}; +int bitpix[] = {8, 16, 32, -32, -64, 0}; + +//char *cmptype[] = {"GZIP_1", NULL}; +//int bitpix[] = {16, 0}; + +static int NX = 10; +static int NY = 10; + +int main (int argc, char **argv) { + + plan_tests (3*3*234); + + diag ("libfits imagecomp.c tests"); + + int i, j; + for (i = 0; cmptype[i]; i++) { + for (j = 0; bitpix[j]; j++) { + if (!strcasecmp (cmptype[i], "RICE_1") && (bitpix[j] < 0)) continue; + if (!strcasecmp (cmptype[i], "RICE_ONE") && (bitpix[j] < 0)) continue; + test_compress (bitpix[j], cmptype[i]); + test_compress_fulltile (bitpix[j], cmptype[i]); + test_compress_fullrange (bitpix[j], cmptype[i]); + } + } + + NX = 100; NY = 33; + for (i = 0; cmptype[i]; i++) { + for (j = 0; bitpix[j]; j++) { + if (!strcasecmp (cmptype[i], "RICE_1") && (bitpix[j] < 0)) continue; + if (!strcasecmp (cmptype[i], "RICE_ONE") && (bitpix[j] < 0)) continue; + test_compress (bitpix[j], cmptype[i]); + test_compress_fulltile (bitpix[j], cmptype[i]); + test_compress_fullrange (bitpix[j], cmptype[i]); + } + } + + NX = 357; NY = 245; + for (i = 0; cmptype[i]; i++) { + for (j = 0; bitpix[j]; j++) { + if (!strcasecmp (cmptype[i], "RICE_1") && (bitpix[j] < 0)) continue; + if (!strcasecmp (cmptype[i], "RICE_ONE") && (bitpix[j] < 0)) continue; + test_compress (bitpix[j], cmptype[i]); + test_compress_fulltile (bitpix[j], cmptype[i]); + test_compress_fullrange (bitpix[j], cmptype[i]); + } + } + return exit_status(); +} + +int test_compress (int bitpix, char *zcmptype) { // make a table, compress, uncompress, compare : use compression zcmptype + + Header rawheader; + Matrix rawmatrix; + Header outheader; + Matrix outmatrix; + + Header theader; + FTable ftable; + ftable.header = &theader; + + diag ("--- starting test_compress with bitpix %d, zcmptype %s ---", bitpix, zcmptype); + ok (gfits_init_header (&rawheader), "init'ed the raw header"); + ok (gfits_init_matrix (&rawmatrix), "init'ed the raw matrix"); + ok (gfits_init_header (&outheader), "init'ed the out header"); + ok (gfits_init_matrix (&outmatrix), "init'ed the out matrix"); + + /* assign the necessary internal values */ + rawheader.bitpix = bitpix; + rawheader.Naxes = 2; + rawheader.Naxis[0] = NX; + rawheader.Naxis[1] = NY; + + /* create the appropriate header and matrix */ + ok (gfits_create_header (&rawheader), "created header"); + ok (gfits_create_matrix (&rawheader, &rawmatrix), "created matrix"); + + char *rawdata_char = (char *) rawmatrix.buffer; + short *rawdata_short = (short *) rawmatrix.buffer; + int *rawdata_int = (int *) rawmatrix.buffer; + float *rawdata_float = (float *) rawmatrix.buffer; + double *rawdata_double = (double *) rawmatrix.buffer; + + int ix, iy; + for (ix = 0; ix < NX; ix++) { + for (iy = 0; iy < NY; iy++) { + int IY = iy + 1; + switch (bitpix) { + case 8: rawdata_char [ix + NX*iy] = ix + IY*IY; break; + case 16: rawdata_short [ix + NX*iy] = ix + IY*IY; break; + case 32: rawdata_int [ix + NX*iy] = ix + IY*IY; break; + case -32: rawdata_float [ix + NX*iy] = ix + IY*IY; break; + case -64: rawdata_double[ix + NX*iy] = ix + IY*IY; break; + default: myAbort ("bad bitpix value"); + } + } + } + + // int Ztile[2] = {NX, NY}; + // ok (gfits_compress_image (&rawheader, &rawmatrix, &ftable, (int *) &Ztile, zcmptype), "compressed image"); + ok (gfits_compress_image (&rawheader, &rawmatrix, &ftable, NULL, zcmptype), "compressed image"); + ok (gfits_uncompress_image (&outheader, &outmatrix, &ftable), "uncompressed image"); + + gfits_free_header (&theader); + gfits_free_table (&ftable); + + int Nbad = 0; + + char *outdata_char = (char *) outmatrix.buffer; + short *outdata_short = (short *) outmatrix.buffer; + int *outdata_int = (int *) outmatrix.buffer; + float *outdata_float = (float *) outmatrix.buffer; + double *outdata_double = (double *) outmatrix.buffer; + + for (ix = 0; ix < NX; ix++) { + for (iy = 0; iy < NY; iy++) { + switch (bitpix) { + case 8: if (rawdata_char [ix + NX*iy] != outdata_char [ix + NX*iy]) Nbad ++; break; + case 16: if (rawdata_short [ix + NX*iy] != outdata_short [ix + NX*iy]) Nbad ++; break; + case 32: if (rawdata_int [ix + NX*iy] != outdata_int [ix + NX*iy]) Nbad ++; break; + case -32: if (rawdata_float [ix + NX*iy] != outdata_float [ix + NX*iy]) Nbad ++; break; + case -64: if (rawdata_double[ix + NX*iy] != outdata_double[ix + NX*iy]) Nbad ++; break; + default: myAbort ("bad bitpix value"); + } + } + } + + ok (!Nbad, "all image pixels match"); + + gfits_free_header (&rawheader); + gfits_free_matrix (&rawmatrix); + + gfits_free_header (&outheader); + gfits_free_matrix (&outmatrix); + return TRUE; +} + +int test_compress_fulltile (int bitpix, char *zcmptype) { // make a table, compress, uncompress, compare : use compression zcmptype + + Header rawheader; + Matrix rawmatrix; + Header outheader; + Matrix outmatrix; + + Header theader; + FTable ftable; + ftable.header = &theader; + + diag ("--- starting test_compress with bitpix %d, zcmptype %s ---", bitpix, zcmptype); + ok (gfits_init_header (&rawheader), "init'ed the raw header"); + ok (gfits_init_matrix (&rawmatrix), "init'ed the raw matrix"); + ok (gfits_init_header (&outheader), "init'ed the out header"); + ok (gfits_init_matrix (&outmatrix), "init'ed the out matrix"); + + /* assign the necessary internal values */ + rawheader.bitpix = bitpix; + rawheader.Naxes = 2; + rawheader.Naxis[0] = NX; + rawheader.Naxis[1] = NY; + + /* create the appropriate header and matrix */ + ok (gfits_create_header (&rawheader), "created header"); + ok (gfits_create_matrix (&rawheader, &rawmatrix), "created matrix"); + + char *rawdata_char = (char *) rawmatrix.buffer; + short *rawdata_short = (short *) rawmatrix.buffer; + int *rawdata_int = (int *) rawmatrix.buffer; + float *rawdata_float = (float *) rawmatrix.buffer; + double *rawdata_double = (double *) rawmatrix.buffer; + + int ix, iy; + for (ix = 0; ix < NX; ix++) { + for (iy = 0; iy < NY; iy++) { + int IY = iy + 1; + switch (bitpix) { + case 8: rawdata_char [ix + NX*iy] = ix + IY*IY; break; + case 16: rawdata_short [ix + NX*iy] = ix + IY*IY; break; + case 32: rawdata_int [ix + NX*iy] = ix + IY*IY; break; + case -32: rawdata_float [ix + NX*iy] = ix + IY*IY; break; + case -64: rawdata_double[ix + NX*iy] = ix + IY*IY; break; + default: myAbort ("bad bitpix value"); + } + } + } + + int Ztile[2] = {NX, NY}; + ok (gfits_compress_image (&rawheader, &rawmatrix, &ftable, (int *) &Ztile, zcmptype), "compressed image"); + // ok (gfits_compress_image (&rawheader, &rawmatrix, &ftable, NULL, zcmptype), "compressed image"); + ok (gfits_uncompress_image (&outheader, &outmatrix, &ftable), "uncompressed image"); + + gfits_free_header (&theader); + gfits_free_table (&ftable); + + int Nbad = 0; + + char *outdata_char = (char *) outmatrix.buffer; + short *outdata_short = (short *) outmatrix.buffer; + int *outdata_int = (int *) outmatrix.buffer; + float *outdata_float = (float *) outmatrix.buffer; + double *outdata_double = (double *) outmatrix.buffer; + + for (ix = 0; ix < NX; ix++) { + for (iy = 0; iy < NY; iy++) { + switch (bitpix) { + case 8: if (rawdata_char [ix + NX*iy] != outdata_char [ix + NX*iy]) Nbad ++; break; + case 16: if (rawdata_short [ix + NX*iy] != outdata_short [ix + NX*iy]) Nbad ++; break; + case 32: if (rawdata_int [ix + NX*iy] != outdata_int [ix + NX*iy]) Nbad ++; break; + case -32: if (rawdata_float [ix + NX*iy] != outdata_float [ix + NX*iy]) Nbad ++; break; + case -64: if (rawdata_double[ix + NX*iy] != outdata_double[ix + NX*iy]) Nbad ++; break; + default: myAbort ("bad bitpix value"); + } + } + } + + ok (!Nbad, "all image pixels match"); + + gfits_free_header (&rawheader); + gfits_free_matrix (&rawmatrix); + + gfits_free_header (&outheader); + gfits_free_matrix (&outmatrix); + return TRUE; +} + +int test_compress_fullrange (int bitpix, char *zcmptype) { // make a table, compress, uncompress, compare : use compression zcmptype + + Header rawheader; + Matrix rawmatrix; + Header outheader; + Matrix outmatrix; + + Header theader; + FTable ftable; + ftable.header = &theader; + + diag ("--- starting test_compress with bitpix %d, zcmptype %s ---", bitpix, zcmptype); + ok (gfits_init_header (&rawheader), "init'ed the raw header"); + ok (gfits_init_matrix (&rawmatrix), "init'ed the raw matrix"); + ok (gfits_init_header (&outheader), "init'ed the out header"); + ok (gfits_init_matrix (&outmatrix), "init'ed the out matrix"); + + /* assign the necessary internal values */ + rawheader.bitpix = bitpix; + rawheader.Naxes = 2; + rawheader.Naxis[0] = NX; + rawheader.Naxis[1] = NY; + + /* create the appropriate header and matrix */ + ok (gfits_create_header (&rawheader), "created header"); + ok (gfits_create_matrix (&rawheader, &rawmatrix), "created matrix"); + + char *rawdata_char = (char *) rawmatrix.buffer; + short *rawdata_short = (short *) rawmatrix.buffer; + int *rawdata_int = (int *) rawmatrix.buffer; + float *rawdata_float = (float *) rawmatrix.buffer; + double *rawdata_double = (double *) rawmatrix.buffer; + + long A = time(NULL); + srand48(A); + // srand48(1); + + int ix, iy; + for (ix = 0; ix < NX; ix++) { + for (iy = 0; iy < NY; iy++) { + switch (bitpix) { + case 8: rawdata_char [ix + NX*iy] = 0xff & lrand48(); break; + case 16: rawdata_short [ix + NX*iy] = 0xffff & lrand48(); break; + case 32: rawdata_int [ix + NX*iy] = 0xffffffff & lrand48(); break; + case -32: rawdata_float [ix + NX*iy] = FLT_MAX * (2.0*drand48() - 1.0); break; + case -64: rawdata_double[ix + NX*iy] = DBL_MAX * (2.0*drand48() - 1.0); break; + default: myAbort ("bad bitpix value"); + } + } + } + + // int Ztile[2] = {NX, NY}; + ok (gfits_compress_image (&rawheader, &rawmatrix, &ftable, NULL, zcmptype), "compressed image"); + + ok (gfits_uncompress_image (&outheader, &outmatrix, &ftable), "uncompressed image"); + + gfits_free_header (&theader); + gfits_free_table (&ftable); + + int Nbad = 0; + + char *outdata_char = (char *) outmatrix.buffer; + short *outdata_short = (short *) outmatrix.buffer; + int *outdata_int = (int *) outmatrix.buffer; + float *outdata_float = (float *) outmatrix.buffer; + double *outdata_double = (double *) outmatrix.buffer; + + for (ix = 0; ix < NX; ix++) { + for (iy = 0; iy < NY; iy++) { + switch (bitpix) { + case 8: if (rawdata_char [ix + NX*iy] != outdata_char [ix + NX*iy]) Nbad ++; break; + case 16: if (rawdata_short [ix + NX*iy] != outdata_short [ix + NX*iy]) Nbad ++; break; + case 32: if (rawdata_int [ix + NX*iy] != outdata_int [ix + NX*iy]) Nbad ++; break; + case -32: if (rawdata_float [ix + NX*iy] != outdata_float [ix + NX*iy]) Nbad ++; break; + case -64: if (rawdata_double[ix + NX*iy] != outdata_double[ix + NX*iy]) Nbad ++; break; + default: myAbort ("bad bitpix value"); + } + } + } + + ok (!Nbad, "all image pixels match"); + + gfits_free_header (&rawheader); + gfits_free_matrix (&rawmatrix); + + gfits_free_header (&outheader); + gfits_free_matrix (&outmatrix); + return TRUE; +} + + +// fprintf (stderr, "%2d %2d : 0x%08x vs 0x%08x vs 0x%08x : %d\n", ix, iy, rawdata_int [ix + NX*iy], outdata_int [ix + NX*iy], tmpbuff[ix + NX*iy], outdata_int [ix + NX*iy] - rawdata_int [ix + NX*iy]); Index: trunk/Ohana/src/libfits/test/ricetest.c =================================================================== --- trunk/Ohana/src/libfits/test/ricetest.c (revision 38441) +++ trunk/Ohana/src/libfits/test/ricetest.c (revision 38441) @@ -0,0 +1,98 @@ +# include +# include +# include "tap_ohana.h" + +# define NPIX 100 +# define NBYTE 4*NPIX + +int main (int argc, char **argv) { + + plan_tests (45); + + diag ("libfits ricetest.c tests"); + + // buffers to store max array + char rawdata[NBYTE]; + char cmpdata[NBYTE]; + char outdata[NBYTE]; + + // ok (1, "failure"); + + if (1) { + char *rawvalue = (char *) rawdata; + char *outvalue = (char *) outdata; + + int j; + for (j = 0; j < 5; j++) { + int i; + for (i = 0; i < NPIX; i++) { + rawvalue[i] = i; + } + + int Ncmp = fits_rcomp_byte (rawvalue, NPIX, cmpdata, NBYTE, 32); + ok (Ncmp > 0, "compressed byte data"); + + int status = fits_rdecomp_byte (cmpdata, Ncmp, outdata, NPIX, 32); + ok (!status, "decompressed byte data"); + + int Nbad = 0; + for (i = 0; i < NPIX; i++) { + if (rawvalue[i] != outvalue[i]) Nbad ++; + } + + ok (!Nbad, "values match"); + } + } + + if (1) { + short *rawvalue = (short *) rawdata; + short *outvalue = (short *) outdata; + + int j; + for (j = 0; j < 5; j++) { + int i; + for (i = 0; i < NPIX; i++) { + rawvalue[i] = i; + } + + int Ncmp = fits_rcomp_short (rawvalue, NPIX, cmpdata, NBYTE, 32); + ok (Ncmp > 0, "compressed short data"); + + int status = fits_rdecomp_short (cmpdata, Ncmp, outdata, NPIX, 32); + ok (!status, "decompressed short data"); + + int Nbad = 0; + for (i = 0; i < NPIX; i++) { + if (rawvalue[i] != outvalue[i]) Nbad ++; + } + + ok (!Nbad, "values match"); + } + } + + if (1) { + int *rawvalue = (int *) rawdata; + int *outvalue = (int *) outdata; + + int j; + for (j = 0; j < 5; j++) { + int i; + for (i = 0; i < NPIX; i++) { + rawvalue[i] = i; + } + + int Ncmp = fits_rcomp (rawvalue, NPIX, cmpdata, NBYTE, 32); + ok (Ncmp > 0, "compressed int data"); + + int status = fits_rdecomp (cmpdata, Ncmp, outdata, NPIX, 32); + ok (!status, "decompressed int data"); + + int Nbad = 0; + for (i = 0; i < NPIX; i++) { + if (rawvalue[i] != outvalue[i]) Nbad ++; + } + ok (!Nbad, "values match"); + } + } + return exit_status(); +} Index: trunk/Ohana/src/libfits/test/tablecomp.c =================================================================== --- trunk/Ohana/src/libfits/test/tablecomp.c (revision 38441) +++ trunk/Ohana/src/libfits/test/tablecomp.c (revision 38441) @@ -0,0 +1,453 @@ +# include +# include +# include "tap_ohana.h" + +int test_compress (char *zcmptype); +int test_compress_empty (char *zcmptype); +int test_compress_fullrange (char *zcmptype); +int test_compress_single_full (char *zcmptype); + +char *cmptype[] = {"NONE", "NONE_1", "NONE_2", "GZIP_1", "GZIP_2", "RICE_1", "RICE_ONE", "AUTO", NULL}; + +int main (int argc, char **argv) { + + plan_tests (668); + + diag ("libfits tablecomp.c tests"); + + // test_compress_empty ("GZIP_1"); + // exit (0); + + test_compress (NULL); + test_compress_fullrange (NULL); + + int i; + for (i = 0; cmptype[i]; i++) { + test_compress (cmptype[i]); + test_compress_fullrange (cmptype[i]); + test_compress_empty (cmptype[i]); + ok (ohana_memcheck_func (TRUE), "no memory corruption"); + } + + return exit_status(); +} + +int test_compress (char *zcmptype) { // test 2: make a table, compress, uncompress, compare : use compression "NONE" + + Header header; + FTable ftable; + + diag ("--- starting test_compress with zcmptype %s ---", zcmptype); + ok (gfits_init_header (&header), "inited the header"); + ok (gfits_init_table (&ftable), "inited the table"); + + ok (gfits_create_table_header (&header, "BINTABLE", "TESTDATA"), "created the table header"); + + ok (gfits_define_bintable_column (&header, "B", "VAL_B", "byte", "none", 1.0, 0.0), "defined byte column"); + ok (gfits_define_bintable_column (&header, "I", "VAL_I", "short", "none", 1.0, 0.0), "defined short column"); + ok (gfits_define_bintable_column (&header, "J", "VAL_J", "int", "none", 1.0, 0.0), "defined int column"); + ok (gfits_define_bintable_column (&header, "K", "VAL_K", "long", "none", 1.0, 0.0), "defined long column"); + ok (gfits_define_bintable_column (&header, "E", "VAL_E", "float", "degree", 1.0, 0.0), "defined float column"); + ok (gfits_define_bintable_column (&header, "D", "VAL_D", "double", "degree", 1.0, 0.0), "defined double column"); + + // generate the output array that carries the data + ok (gfits_create_table (&header, &ftable), "created the basic table"); + + int Nval = 1000; + char *VAL_B; ALLOCATE (VAL_B, char, Nval); + short *VAL_I; ALLOCATE (VAL_I, short, Nval); + int *VAL_J; ALLOCATE (VAL_J, int, Nval); + int64_t *VAL_K; ALLOCATE (VAL_K, int64_t, Nval); + float *VAL_E; ALLOCATE (VAL_E, float, Nval); + double *VAL_D; ALLOCATE (VAL_D, double, Nval); + + int i; + for (i = 0; i < Nval; i++) { + VAL_B[i] = (i % 255) - 128; + VAL_I[i] = i; + VAL_J[i] = 10000*i + 100*i; + VAL_K[i] = 10000*i + 330*i; + VAL_E[i] = 0.1*i; + VAL_D[i] = 1.001*i; + } + + // add the columns to the output array + ok (gfits_set_bintable_column (&header, &ftable, "VAL_B", VAL_B, Nval), "set byte table column"); + ok (gfits_set_bintable_column (&header, &ftable, "VAL_I", VAL_I, Nval), "set short table column"); + ok (gfits_set_bintable_column (&header, &ftable, "VAL_J", VAL_J, Nval), "set int table column"); + ok (gfits_set_bintable_column (&header, &ftable, "VAL_K", VAL_K, Nval), "set long table column"); + ok (gfits_set_bintable_column (&header, &ftable, "VAL_E", VAL_E, Nval), "set float table column"); + ok (gfits_set_bintable_column (&header, &ftable, "VAL_D", VAL_D, Nval), "set double table column"); + + Header *outheader = &header; + FTable *outtable = &ftable; + + if (zcmptype) { + Header cmpheader; + FTable cmptable; + cmptable.header = &cmpheader; + ok (gfits_compress_table (&ftable, &cmptable, 10, zcmptype), "compressed table"); + + Header rawheader; + FTable rawtable; + rawtable.header = &rawheader; + ok (gfits_uncompress_table (&cmptable, &rawtable), "uncompressed table"); + + gfits_free_header (&cmpheader); + gfits_free_table (&cmptable); + + outheader = &rawheader; + outtable = &rawtable; + } + + char type[16]; + int Ncol; + off_t Nrow; + + char *VAL_B_t = gfits_get_bintable_column_data (outheader, outtable, "VAL_B", type, &Nrow, &Ncol); ok (!strcmp (type, "byte"), "read byte table column"); ok (Nrow == Nval, "right number of rows"); ok (Ncol == 1, "right number of cols"); + short *VAL_I_t = gfits_get_bintable_column_data (outheader, outtable, "VAL_I", type, &Nrow, &Ncol); ok (!strcmp (type, "short"), "read short table column"); ok (Nrow == Nval, "right number of rows"); ok (Ncol == 1, "right number of cols"); + int *VAL_J_t = gfits_get_bintable_column_data (outheader, outtable, "VAL_J", type, &Nrow, &Ncol); ok (!strcmp (type, "int"), "read int table column"); ok (Nrow == Nval, "right number of rows"); ok (Ncol == 1, "right number of cols"); + int64_t *VAL_K_t = gfits_get_bintable_column_data (outheader, outtable, "VAL_K", type, &Nrow, &Ncol); ok (!strcmp (type, "int64_t"), "read int64_t table column"); ok (Nrow == Nval, "right number of rows"); ok (Ncol == 1, "right number of cols"); + float *VAL_E_t = gfits_get_bintable_column_data (outheader, outtable, "VAL_E", type, &Nrow, &Ncol); ok (!strcmp (type, "float"), "read float table column"); ok (Nrow == Nval, "right number of rows"); ok (Ncol == 1, "right number of cols"); + double *VAL_D_t = gfits_get_bintable_column_data (outheader, outtable, "VAL_D", type, &Nrow, &Ncol); ok (!strcmp (type, "double"), "read double table column"); ok (Nrow == Nval, "right number of rows"); ok (Ncol == 1, "right number of cols"); + + // count mismatched values + int NVAL_B_bad = 0; + int NVAL_I_bad = 0; + int NVAL_J_bad = 0; + int NVAL_K_bad = 0; + int NVAL_E_bad = 0; + int NVAL_D_bad = 0; + + // for (i = 0; i < Nval; i++) { + // fprintf (stderr, "%0llx : %0llx\n", (long long) VAL_K[i], (long long) VAL_K_t[i]); + // } + + for (i = 0; i < Nval; i++) { + if (VAL_B[i] != VAL_B_t[i]) NVAL_B_bad++; + if (VAL_I[i] != VAL_I_t[i]) NVAL_I_bad++; + if (VAL_J[i] != VAL_J_t[i]) NVAL_J_bad++; + if (VAL_K[i] != VAL_K_t[i]) NVAL_K_bad++; + if (VAL_E[i] != VAL_E_t[i]) NVAL_E_bad++; + if (VAL_D[i] != VAL_D_t[i]) NVAL_D_bad++; + } + + ok (!NVAL_B_bad, "byte values match (input vs output)"); + ok (!NVAL_I_bad, "short values match (input vs output)"); + ok (!NVAL_J_bad, "int values match (input vs output)"); + ok (!NVAL_K_bad, "int64_t values match (input vs output)"); + ok (!NVAL_E_bad, "float values match (input vs output)"); + ok (!NVAL_D_bad, "double values match (input vs output)"); + + gfits_free_header (&header); + gfits_free_table (&ftable); + if (zcmptype) { + gfits_free_header (outheader); + gfits_free_table (outtable); + } + return TRUE; +} + +int test_compress_fullrange (char *zcmptype) { // test 2: make a table, compress, uncompress, compare : use compression "NONE" + + Header header; + FTable ftable; + + diag ("--- starting test_compress with zcmptype %s ---", zcmptype); + ok (gfits_init_header (&header), "inited the header"); + ok (gfits_init_table (&ftable), "inited the table"); + + ok (gfits_create_table_header (&header, "BINTABLE", "TESTDATA"), "created the table header"); + + ok (gfits_define_bintable_column (&header, "B", "VAL_B", "byte", "none", 1.0, 0.0), "defined byte column"); + ok (gfits_define_bintable_column (&header, "I", "VAL_I", "short", "none", 1.0, 0.0), "defined short column"); + ok (gfits_define_bintable_column (&header, "J", "VAL_J", "int", "none", 1.0, 0.0), "defined int column"); + ok (gfits_define_bintable_column (&header, "K", "VAL_K", "long", "none", 1.0, 0.0), "defined long column"); + ok (gfits_define_bintable_column (&header, "E", "VAL_E", "float", "degree", 1.0, 0.0), "defined float column"); + ok (gfits_define_bintable_column (&header, "D", "VAL_D", "double", "degree", 1.0, 0.0), "defined double column"); + + // generate the output array that carries the data + ok (gfits_create_table (&header, &ftable), "created the basic table"); + + int Nval = 1000; + char *VAL_B; ALLOCATE (VAL_B, char, Nval); + short *VAL_I; ALLOCATE (VAL_I, short, Nval); + int *VAL_J; ALLOCATE (VAL_J, int, Nval); + int64_t *VAL_K; ALLOCATE (VAL_K, int64_t, Nval); + float *VAL_E; ALLOCATE (VAL_E, float, Nval); + double *VAL_D; ALLOCATE (VAL_D, double, Nval); + + long A = time(NULL); + srand48(A); + + int i; + for (i = 0; i < Nval; i++) { + VAL_B[i] = 0xff & lrand48(); + VAL_I[i] = 0xffff & lrand48(); + VAL_J[i] = 0xffffffff & lrand48(); + VAL_K[i] = (((int64_t)lrand48()) << 32) + (int64_t)lrand48(); + VAL_E[i] = FLT_MAX * (2.0*drand48() - 1.0); + VAL_D[i] = DBL_MAX * (2.0*drand48() - 1.0); + } + + // add the columns to the output array + ok (gfits_set_bintable_column (&header, &ftable, "VAL_B", VAL_B, Nval), "set byte table column"); + ok (gfits_set_bintable_column (&header, &ftable, "VAL_I", VAL_I, Nval), "set short table column"); + ok (gfits_set_bintable_column (&header, &ftable, "VAL_J", VAL_J, Nval), "set int table column"); + ok (gfits_set_bintable_column (&header, &ftable, "VAL_K", VAL_K, Nval), "set long table column"); + ok (gfits_set_bintable_column (&header, &ftable, "VAL_E", VAL_E, Nval), "set float table column"); + ok (gfits_set_bintable_column (&header, &ftable, "VAL_D", VAL_D, Nval), "set double table column"); + + Header *outheader = &header; + FTable *outtable = &ftable; + + if (zcmptype) { + Header cmpheader; + FTable cmptable; + cmptable.header = &cmpheader; + ok (gfits_compress_table (&ftable, &cmptable, 10, zcmptype), "compressed table"); + + Header rawheader; + FTable rawtable; + rawtable.header = &rawheader; + ok (gfits_uncompress_table (&cmptable, &rawtable), "uncompressed table"); + + gfits_free_header (&cmpheader); + gfits_free_table (&cmptable); + + outheader = &rawheader; + outtable = &rawtable; + } + + char type[16]; + int Ncol; + off_t Nrow; + + char *VAL_B_t = gfits_get_bintable_column_data (outheader, outtable, "VAL_B", type, &Nrow, &Ncol); ok (!strcmp (type, "byte"), "read byte table column"); ok (Nrow == Nval, "right number of rows"); ok (Ncol == 1, "right number of cols"); + short *VAL_I_t = gfits_get_bintable_column_data (outheader, outtable, "VAL_I", type, &Nrow, &Ncol); ok (!strcmp (type, "short"), "read short table column"); ok (Nrow == Nval, "right number of rows"); ok (Ncol == 1, "right number of cols"); + int *VAL_J_t = gfits_get_bintable_column_data (outheader, outtable, "VAL_J", type, &Nrow, &Ncol); ok (!strcmp (type, "int"), "read int table column"); ok (Nrow == Nval, "right number of rows"); ok (Ncol == 1, "right number of cols"); + int64_t *VAL_K_t = gfits_get_bintable_column_data (outheader, outtable, "VAL_K", type, &Nrow, &Ncol); ok (!strcmp (type, "int64_t"), "read int64_t table column"); ok (Nrow == Nval, "right number of rows"); ok (Ncol == 1, "right number of cols"); + float *VAL_E_t = gfits_get_bintable_column_data (outheader, outtable, "VAL_E", type, &Nrow, &Ncol); ok (!strcmp (type, "float"), "read float table column"); ok (Nrow == Nval, "right number of rows"); ok (Ncol == 1, "right number of cols"); + double *VAL_D_t = gfits_get_bintable_column_data (outheader, outtable, "VAL_D", type, &Nrow, &Ncol); ok (!strcmp (type, "double"), "read double table column"); ok (Nrow == Nval, "right number of rows"); ok (Ncol == 1, "right number of cols"); + + // count mismatched values + int NVAL_B_bad = 0; + int NVAL_I_bad = 0; + int NVAL_J_bad = 0; + int NVAL_K_bad = 0; + int NVAL_E_bad = 0; + int NVAL_D_bad = 0; + + for (i = 0; i < Nval; i++) { + if (VAL_B[i] != VAL_B_t[i]) NVAL_B_bad++; + if (VAL_I[i] != VAL_I_t[i]) NVAL_I_bad++; + if (VAL_J[i] != VAL_J_t[i]) NVAL_J_bad++; + if (VAL_K[i] != VAL_K_t[i]) NVAL_K_bad++; + if (VAL_E[i] != VAL_E_t[i]) NVAL_E_bad++; + if (VAL_D[i] != VAL_D_t[i]) NVAL_D_bad++; + } + + ok (!NVAL_B_bad, "byte values match (input vs output)"); + ok (!NVAL_I_bad, "short values match (input vs output)"); + ok (!NVAL_J_bad, "int values match (input vs output)"); + ok (!NVAL_K_bad, "int64_t values match (input vs output)"); + ok (!NVAL_E_bad, "float values match (input vs output)"); + ok (!NVAL_D_bad, "double values match (input vs output)"); + + gfits_free_header (&header); + gfits_free_table (&ftable); + if (zcmptype) { + gfits_free_header (outheader); + gfits_free_table (outtable); + } + return TRUE; +} + +int test_compress_single_full (char *zcmptype) { // test 2: make a table, compress, uncompress, compare : use compression "NONE" + + Header header; + FTable ftable; + + diag ("--- starting test_compress with zcmptype %s ---", zcmptype); + ok (gfits_init_header (&header), "inited the header"); + ok (gfits_init_table (&ftable), "inited the table"); + + ok (gfits_create_table_header (&header, "BINTABLE", "TESTDATA"), "created the table header"); + + ok (gfits_define_bintable_column (&header, "I", "SEQ", "seq", "none", 1.0, 0.0), "defined short column"); + + // generate the output array that carries the data + ok (gfits_create_table (&header, &ftable), "created the basic table"); + + int Nval = 1000; + short *SEQ; ALLOCATE (SEQ, short, Nval); + + int i; + for (i = 0; i < Nval; i++) { + SEQ[i] = i; + } + + // add the columns to the output array + ok (gfits_set_bintable_column (&header, &ftable, "SEQ", SEQ, Nval), "set short table column"); + + Header *outheader = &header; + FTable *outtable = &ftable; + + if (zcmptype) { + Header cmpheader; + FTable cmptable; + cmptable.header = &cmpheader; + ok (gfits_compress_table (&ftable, &cmptable, 0, zcmptype), "compressed table"); + + Header rawheader; + FTable rawtable; + rawtable.header = &rawheader; + ok (gfits_uncompress_table (&cmptable, &rawtable), "uncompressed table"); + + gfits_free_header (&cmpheader); + gfits_free_table (&cmptable); + + outheader = &rawheader; + outtable = &rawtable; + } + + char type[16]; + int Ncol; + off_t Nrow; + + short *SEQ_t = gfits_get_bintable_column_data (outheader, outtable, "SEQ", type, &Nrow, &Ncol); ok (!strcmp (type, "short"), "read short table column"); ok (Nrow == Nval, "right number of rows"); ok (Ncol == 1, "right number of cols"); + + // count mismatched values + int Nseq = 0; + for (i = 0; i < Nval; i++) { + if (SEQ[i] != SEQ_t[i]) Nseq++; + } + + ok (!Nseq, "short values match (input vs output)"); + + gfits_free_header (&header); + gfits_free_table (&ftable); + if (zcmptype) { + gfits_free_header (outheader); + gfits_free_table (outtable); + } + return TRUE; +} + +int test_compress_empty (char *zcmptype) { // test 2: make a table, compress, uncompress, compare : use compression "NONE" + + Header header; + FTable ftable; + + diag ("--- starting test_compress with zcmptype %s ---", zcmptype); + ok (gfits_init_header (&header), "inited the header"); + ok (gfits_init_table (&ftable), "inited the table"); + + ok (gfits_create_table_header (&header, "BINTABLE", "TESTDATA"), "created the table header"); + + ok (gfits_define_bintable_column (&header, "B", "VAL_B", "byte", "none", 1.0, 0.0), "defined byte column"); + ok (gfits_define_bintable_column (&header, "I", "VAL_I", "short", "none", 1.0, 0.0), "defined short column"); + ok (gfits_define_bintable_column (&header, "J", "VAL_J", "int", "none", 1.0, 0.0), "defined int column"); + ok (gfits_define_bintable_column (&header, "K", "VAL_K", "long", "none", 1.0, 0.0), "defined long column"); + ok (gfits_define_bintable_column (&header, "E", "VAL_E", "float", "degree", 1.0, 0.0), "defined float column"); + ok (gfits_define_bintable_column (&header, "D", "VAL_D", "double", "degree", 1.0, 0.0), "defined double column"); + + // generate the output array that carries the data + ok (gfits_create_table (&header, &ftable), "created the basic table"); + + int Nval = 0; + char *VAL_B; ALLOCATE (VAL_B, char, Nval); + short *VAL_I; ALLOCATE (VAL_I, short, Nval); + int *VAL_J; ALLOCATE (VAL_J, int, Nval); + int64_t *VAL_K; ALLOCATE (VAL_K, int64_t, Nval); + float *VAL_E; ALLOCATE (VAL_E, float, Nval); + double *VAL_D; ALLOCATE (VAL_D, double, Nval); + + int i; + for (i = 0; i < Nval; i++) { + VAL_B[i] = (i % 255) - 128; + VAL_I[i] = i; + VAL_J[i] = 10000*i + 100*i; + VAL_K[i] = 10000*i + 330*i; + VAL_E[i] = 0.1*i; + VAL_D[i] = 1.001*i; + } + + // add the columns to the output array + ok (gfits_set_bintable_column (&header, &ftable, "VAL_B", VAL_B, Nval), "set byte table column"); + ok (gfits_set_bintable_column (&header, &ftable, "VAL_I", VAL_I, Nval), "set short table column"); + ok (gfits_set_bintable_column (&header, &ftable, "VAL_J", VAL_J, Nval), "set int table column"); + ok (gfits_set_bintable_column (&header, &ftable, "VAL_K", VAL_K, Nval), "set long table column"); + ok (gfits_set_bintable_column (&header, &ftable, "VAL_E", VAL_E, Nval), "set float table column"); + ok (gfits_set_bintable_column (&header, &ftable, "VAL_D", VAL_D, Nval), "set double table column"); + + Header *outheader = &header; + FTable *outtable = &ftable; + + if (zcmptype) { + Header cmpheader; + FTable cmptable; + cmptable.header = &cmpheader; + ok (gfits_compress_table (&ftable, &cmptable, 0, zcmptype), "compressed table"); + + Header rawheader; + FTable rawtable; + rawtable.header = &rawheader; + ok (gfits_uncompress_table (&cmptable, &rawtable), "uncompressed table"); + + gfits_free_header (&cmpheader); + gfits_free_table (&cmptable); + + outheader = &rawheader; + outtable = &rawtable; + } + + char type[16]; + int Ncol; + off_t Nrow; + + char *VAL_B_t = gfits_get_bintable_column_data (outheader, outtable, "VAL_B", type, &Nrow, &Ncol); ok (!strcmp (type, "byte"), "read byte table column"); ok (Nrow == Nval, "right number of rows"); ok (Ncol == 1, "right number of cols"); + short *VAL_I_t = gfits_get_bintable_column_data (outheader, outtable, "VAL_I", type, &Nrow, &Ncol); ok (!strcmp (type, "short"), "read short table column"); ok (Nrow == Nval, "right number of rows"); ok (Ncol == 1, "right number of cols"); + int *VAL_J_t = gfits_get_bintable_column_data (outheader, outtable, "VAL_J", type, &Nrow, &Ncol); ok (!strcmp (type, "int"), "read int table column"); ok (Nrow == Nval, "right number of rows"); ok (Ncol == 1, "right number of cols"); + int64_t *VAL_K_t = gfits_get_bintable_column_data (outheader, outtable, "VAL_K", type, &Nrow, &Ncol); ok (!strcmp (type, "int64_t"), "read int64_t table column"); ok (Nrow == Nval, "right number of rows"); ok (Ncol == 1, "right number of cols"); + float *VAL_E_t = gfits_get_bintable_column_data (outheader, outtable, "VAL_E", type, &Nrow, &Ncol); ok (!strcmp (type, "float"), "read float table column"); ok (Nrow == Nval, "right number of rows"); ok (Ncol == 1, "right number of cols"); + double *VAL_D_t = gfits_get_bintable_column_data (outheader, outtable, "VAL_D", type, &Nrow, &Ncol); ok (!strcmp (type, "double"), "read double table column"); ok (Nrow == Nval, "right number of rows"); ok (Ncol == 1, "right number of cols"); + + // count mismatched values + int NVAL_B_bad = 0; + int NVAL_I_bad = 0; + int NVAL_J_bad = 0; + int NVAL_K_bad = 0; + int NVAL_E_bad = 0; + int NVAL_D_bad = 0; + + // for (i = 0; i < Nval; i++) { + // fprintf (stderr, "%0llx : %0llx\n", (long long) VAL_K[i], (long long) VAL_K_t[i]); + // } + + for (i = 0; i < Nval; i++) { + if (VAL_B[i] != VAL_B_t[i]) NVAL_B_bad++; + if (VAL_I[i] != VAL_I_t[i]) NVAL_I_bad++; + if (VAL_J[i] != VAL_J_t[i]) NVAL_J_bad++; + if (VAL_K[i] != VAL_K_t[i]) NVAL_K_bad++; + if (VAL_E[i] != VAL_E_t[i]) NVAL_E_bad++; + if (VAL_D[i] != VAL_D_t[i]) NVAL_D_bad++; + } + + ok (!NVAL_B_bad, "byte values match (input vs output)"); + ok (!NVAL_I_bad, "short values match (input vs output)"); + ok (!NVAL_J_bad, "int values match (input vs output)"); + ok (!NVAL_K_bad, "int64_t values match (input vs output)"); + ok (!NVAL_E_bad, "float values match (input vs output)"); + ok (!NVAL_D_bad, "double values match (input vs output)"); + + gfits_free_header (&header); + gfits_free_table (&ftable); + if (zcmptype) { + gfits_free_header (outheader); + gfits_free_table (outtable); + } + return TRUE; +} + + Index: trunk/Ohana/src/libfits/test/tap.dvo =================================================================== --- trunk/Ohana/src/libfits/test/tap.dvo (revision 38441) +++ trunk/Ohana/src/libfits/test/tap.dvo (revision 38441) @@ -0,0 +1,116 @@ +# -*-sh-*- + +if (not($?TAP_BREAK)) set TAP_BREAK = 0 +if (not($?TAP_VERBOSE)) set TAP_VERBOSE = 0 + +macro tapEXEC + if ($0 < 3) + echo "USAGE: tapEXEC (command)" + break + end + + local i cmd result + + $cmd = "" + for i 1 $0 + $cmd = $cmd $$i + end + + if ($TAP_VERBOSE >= 2) echo $cmd + + break -auto off + exec $cmd >& tap.log + $result = $STATUS + break -auto on + + if (not($result)) + exec cat tap.log + echo "not ok : $cmd" + $TAP_NFAIL ++ + $TAP_LAST = 0 + if ($TAP_BREAK) + break + end + else + if ($TAP_VERBOSE >= 1) echo "ok : $cmd" + $TAP_LAST = 1 + end + $TAP_NDONE ++ +end + +macro tapOK + if ($0 != 3) + echo "USAGE: tapOK (condition) (message)" + break + end + + if ($1) + if ($TAP_VERBOSE >= 1) + echo "ok : $2" + end + $TAP_LAST = 1 + else + echo "not ok : $2" + $TAP_NFAIL ++ + $TAP_LAST = 0 + if ($TAP_BREAK) + break + end + end + $TAP_NDONE ++ +end + +macro tapPLAN + if ($0 != 2) + echo "USAGE: tapPLAN (Ntests)" + break + end + + $TAP_NTEST = $1 + $TAP_NFAIL = 0 + $TAP_NSKIP = 0 + $TAP_NDONE = 0 + if (not($?TAP_BREAK)) set TAP_BREAK = 0 +end + +macro tapDIAG + if ($0 != 3) + echo "USAGE: tapDIAG (level) (message)" + break + end + + if ($TAP_VERBOSE >= $1) echo "$2" +end + +macro tapSKIP + if ($0 != 2) + echo "USAGE: tapSKIP (nskip)" + break + end + + $TAP_NSKIP += $1 + $TAP_NDONE += $1 +end + +macro tapDONE + if ($0 != 1) + echo "USAGE: tapDONE" + break + end + + if ($TAP_NDONE != $TAP_NTEST) + echo "planned tests ($TAP_NTEST) not equal to done tests ($TAP_NDONE)" + end + + if ($TAP_NFAIL) + echo "failed $TAP_NFAIL of $TAP_NDONE" + end + + if ($TAP_NSKIP) + echo "skipped $TAP_NSKIP of $TAP_NDONE" + end + + if ($TAP_NFAIL + $TAP_NSKIP == 0) + echo "passed $TAP_NDONE tests" + end +end Index: trunk/Ohana/src/libfits/test/tcomptiming.c =================================================================== --- trunk/Ohana/src/libfits/test/tcomptiming.c (revision 38441) +++ trunk/Ohana/src/libfits/test/tcomptiming.c (revision 38441) @@ -0,0 +1,207 @@ +# include +# include +# include "tap_ohana.h" + +int test_compress_timing (char *zcmptype, int Ntile); +void gfits_uncompress_timing (); + +char *cmptype[] = {"NONE", "NONE_1", "NONE_2", "GZIP_1", "GZIP_2", "RICE_1", "RICE_ONE", "AUTO", NULL}; + +static int Nval = 10000000; + +int main (int argc, char **argv) { + + plan_tests (668); + // plan_skip_all ("skipping"); + + diag ("libfits tablecomp.c tests"); + + // test_compress_empty ("GZIP_1"); + // exit (0); + + float lNtileMax = log10(Nval) + 0.1; + float lNtile = 0.0; + for (lNtile = 1.0; lNtile < lNtileMax; lNtile += 1.0) { + int Ntile = pow (10.0, lNtile); + fprintf (stderr, "--- Ntile = %d ---\n", Ntile); + test_compress_timing ("NONE_2", Ntile); + gfits_compress_timing(); + gfits_uncompress_timing(); + } + ok (ohana_memcheck_func (TRUE), "no memory corruption"); + + return exit_status(); +} + +int test_compress_timing (char *zcmptype, int Ntile) { // test 2: make a table, compress, uncompress, compare : use compression "NONE" + + Header header; + FTable ftable; + + diag ("--- starting test_compress with zcmptype %s ---", zcmptype); + ok (gfits_init_header (&header), "inited the header"); + ok (gfits_init_table (&ftable), "inited the table"); + + ok (gfits_create_table_header (&header, "BINTABLE", "TESTDATA"), "created the table header"); + + ok (gfits_define_bintable_column (&header, "B", "VAL_B", "byte", "none", 1.0, 0.0), "defined byte column"); + ok (gfits_define_bintable_column (&header, "I", "VAL_I", "short", "none", 1.0, 0.0), "defined short column"); + ok (gfits_define_bintable_column (&header, "J", "VAL_J", "int", "none", 1.0, 0.0), "defined int column"); + ok (gfits_define_bintable_column (&header, "K", "VAL_K", "long", "none", 1.0, 0.0), "defined long column"); + ok (gfits_define_bintable_column (&header, "E", "VAL_E", "float", "degree", 1.0, 0.0), "defined float column"); + ok (gfits_define_bintable_column (&header, "D", "VAL_D", "double", "degree", 1.0, 0.0), "defined double column"); + + // generate the output array that carries the data + ok (gfits_create_table (&header, &ftable), "created the basic table"); + + // XXX EAM: + struct timeval startTimer, stopTimer; + float dtime; + gettimeofday (&startTimer, (void *) NULL); + + float timeSum1 = 0.0; + float timeSum2 = 0.0; + float timeSum3 = 0.0; + float timeSum4 = 0.0; + float timeSum5 = 0.0; + float timeSum6 = 0.0; + float timeSum7 = 0.0; + + char *VAL_B; ALLOCATE (VAL_B, char, Nval); + short *VAL_I; ALLOCATE (VAL_I, short, Nval); + int *VAL_J; ALLOCATE (VAL_J, int, Nval); + int64_t *VAL_K; ALLOCATE (VAL_K, int64_t, Nval); + float *VAL_E; ALLOCATE (VAL_E, float, Nval); + double *VAL_D; ALLOCATE (VAL_D, double, Nval); + + long A = time(NULL); + srand48(A); + + // XXX TIMER 1 + gettimeofday (&stopTimer, (void *) NULL); + dtime = DTIME (stopTimer, startTimer); + timeSum1 += dtime; + gettimeofday (&startTimer, (void *) NULL); + + int i; + for (i = 0; i < Nval; i++) { + VAL_B[i] = 0xff & lrand48(); + VAL_I[i] = 0xffff & lrand48(); + VAL_J[i] = 0xffffffff & lrand48(); + VAL_K[i] = (((int64_t)lrand48()) << 32) + (int64_t)lrand48(); + VAL_E[i] = FLT_MAX * (2.0*drand48() - 1.0); + VAL_D[i] = DBL_MAX * (2.0*drand48() - 1.0); + } + + // XXX TIMER 1 + gettimeofday (&stopTimer, (void *) NULL); + dtime = DTIME (stopTimer, startTimer); + timeSum2 += dtime; + gettimeofday (&startTimer, (void *) NULL); + + // add the columns to the output array + ok (gfits_set_bintable_column (&header, &ftable, "VAL_B", VAL_B, Nval), "set byte table column"); + ok (gfits_set_bintable_column (&header, &ftable, "VAL_I", VAL_I, Nval), "set short table column"); + ok (gfits_set_bintable_column (&header, &ftable, "VAL_J", VAL_J, Nval), "set int table column"); + ok (gfits_set_bintable_column (&header, &ftable, "VAL_K", VAL_K, Nval), "set long table column"); + ok (gfits_set_bintable_column (&header, &ftable, "VAL_E", VAL_E, Nval), "set float table column"); + ok (gfits_set_bintable_column (&header, &ftable, "VAL_D", VAL_D, Nval), "set double table column"); + + // XXX TIMER 1 + gettimeofday (&stopTimer, (void *) NULL); + dtime = DTIME (stopTimer, startTimer); + timeSum3 += dtime; + gettimeofday (&startTimer, (void *) NULL); + + Header *outheader = &header; + FTable *outtable = &ftable; + + if (zcmptype) { + Header cmpheader; + FTable cmptable; + cmptable.header = &cmpheader; + ok (gfits_compress_table (&ftable, &cmptable, Ntile, zcmptype), "compressed table"); + + Header rawheader; + FTable rawtable; + rawtable.header = &rawheader; + ok (gfits_uncompress_table (&cmptable, &rawtable), "uncompressed table"); + + gfits_free_header (&cmpheader); + gfits_free_table (&cmptable); + + outheader = &rawheader; + outtable = &rawtable; + } + + char type[16]; + int Ncol; + off_t Nrow; + + // XXX TIMER 1 + gettimeofday (&stopTimer, (void *) NULL); + dtime = DTIME (stopTimer, startTimer); + timeSum4 += dtime; + gettimeofday (&startTimer, (void *) NULL); + + char *VAL_B_t = gfits_get_bintable_column_data (outheader, outtable, "VAL_B", type, &Nrow, &Ncol); ok (!strcmp (type, "byte"), "read byte table column"); ok (Nrow == Nval, "right number of rows"); ok (Ncol == 1, "right number of cols"); + short *VAL_I_t = gfits_get_bintable_column_data (outheader, outtable, "VAL_I", type, &Nrow, &Ncol); ok (!strcmp (type, "short"), "read short table column"); ok (Nrow == Nval, "right number of rows"); ok (Ncol == 1, "right number of cols"); + int *VAL_J_t = gfits_get_bintable_column_data (outheader, outtable, "VAL_J", type, &Nrow, &Ncol); ok (!strcmp (type, "int"), "read int table column"); ok (Nrow == Nval, "right number of rows"); ok (Ncol == 1, "right number of cols"); + int64_t *VAL_K_t = gfits_get_bintable_column_data (outheader, outtable, "VAL_K", type, &Nrow, &Ncol); ok (!strcmp (type, "int64_t"), "read int64_t table column"); ok (Nrow == Nval, "right number of rows"); ok (Ncol == 1, "right number of cols"); + float *VAL_E_t = gfits_get_bintable_column_data (outheader, outtable, "VAL_E", type, &Nrow, &Ncol); ok (!strcmp (type, "float"), "read float table column"); ok (Nrow == Nval, "right number of rows"); ok (Ncol == 1, "right number of cols"); + double *VAL_D_t = gfits_get_bintable_column_data (outheader, outtable, "VAL_D", type, &Nrow, &Ncol); ok (!strcmp (type, "double"), "read double table column"); ok (Nrow == Nval, "right number of rows"); ok (Ncol == 1, "right number of cols"); + + // XXX TIMER 1 + gettimeofday (&stopTimer, (void *) NULL); + dtime = DTIME (stopTimer, startTimer); + timeSum5 += dtime; + gettimeofday (&startTimer, (void *) NULL); + + // count mismatched values + int NVAL_B_bad = 0; + int NVAL_I_bad = 0; + int NVAL_J_bad = 0; + int NVAL_K_bad = 0; + int NVAL_E_bad = 0; + int NVAL_D_bad = 0; + + for (i = 0; i < Nval; i++) { + if (VAL_B[i] != VAL_B_t[i]) NVAL_B_bad++; + if (VAL_I[i] != VAL_I_t[i]) NVAL_I_bad++; + if (VAL_J[i] != VAL_J_t[i]) NVAL_J_bad++; + if (VAL_K[i] != VAL_K_t[i]) NVAL_K_bad++; + if (VAL_E[i] != VAL_E_t[i]) NVAL_E_bad++; + if (VAL_D[i] != VAL_D_t[i]) NVAL_D_bad++; + } + + // XXX TIMER 1 + gettimeofday (&stopTimer, (void *) NULL); + dtime = DTIME (stopTimer, startTimer); + timeSum6 += dtime; + gettimeofday (&startTimer, (void *) NULL); + + ok (!NVAL_B_bad, "byte values match (input vs output)"); + ok (!NVAL_I_bad, "short values match (input vs output)"); + ok (!NVAL_J_bad, "int values match (input vs output)"); + ok (!NVAL_K_bad, "int64_t values match (input vs output)"); + ok (!NVAL_E_bad, "float values match (input vs output)"); + ok (!NVAL_D_bad, "double values match (input vs output)"); + + gfits_free_header (&header); + gfits_free_table (&ftable); + if (zcmptype) { + gfits_free_header (outheader); + gfits_free_table (outtable); + } + // XXX TIMER 1 + gettimeofday (&stopTimer, (void *) NULL); + dtime = DTIME (stopTimer, startTimer); + timeSum7 += dtime; + gettimeofday (&startTimer, (void *) NULL); + + fprintf (stderr, "out times: %f %f %f %f %f %f %f\n", timeSum1, timeSum2, timeSum3, timeSum4, timeSum5, timeSum6, timeSum7); + + return TRUE; +} + + Index: trunk/Ohana/src/libfits/test/test.zlib.c =================================================================== --- trunk/Ohana/src/libfits/test/test.zlib.c (revision 38441) +++ trunk/Ohana/src/libfits/test/test.zlib.c (revision 38441) @@ -0,0 +1,73 @@ +# include +# include +# include +# include +# include + +# define NSRC 100 +# define NTGT 200 + +int main (int argc, char **argv) { + + int i, err; + char srcbuf[NSRC], tgtbuf[NTGT], outbuf[NTGT]; + + { + z_stream zdn; + zdn.next_in = srcbuf; + zdn.avail_in = NSRC; + zdn.next_out = tgtbuf; + zdn.avail_out = NTGT; + + zdn.zalloc = Z_NULL; + zdn.zfree = Z_NULL; + zdn.opaque = Z_NULL; + + for (i = 0; i < NSRC; i++) srcbuf[i] = i; + memset (tgtbuf, 0, NTGT); + memset (outbuf, 0, NTGT); + + // the '5' is the compression level: make this a user argument + err = deflateInit(&zdn, 5); + if (err != Z_OK) { fprintf (stderr, "error 1: %d vs %d\n", err, Z_OK); exit (1); } + + err = deflate(&zdn, Z_FINISH); + if (err != Z_STREAM_END) { fprintf (stderr, "error 2: %d vs %d\n", err, Z_STREAM_END); exit (2); } + + assert (zdn.total_out <= NTGT); + + fprintf (stderr, "result: %d bytes\n", (int) zdn.total_out); + for (i = 0; i < zdn.total_out; i++) { + fprintf (stderr, "%d : 0x%02hhx\n", i, tgtbuf[i]); + } + } + + { + z_stream zup; + zup.next_in = tgtbuf; + zup.avail_in = NTGT; + zup.next_out = outbuf; + zup.avail_out = NTGT; + + zup.zalloc = Z_NULL; + zup.zfree = Z_NULL; + zup.opaque = Z_NULL; + + // the '5' is the compression level: make this a user argument + err = inflateInit(&zup); + if (err != Z_OK) { fprintf (stderr, "error 1: %d vs %d\n", err, Z_OK); exit (1); } + + err = inflate(&zup, Z_FINISH); + if (err != Z_STREAM_END) { fprintf (stderr, "error 2: %d vs %d\n", err, Z_STREAM_END); exit (2); } + + assert (zup.total_out <= NTGT); + + fprintf (stderr, "result: %d bytes\n", (int) zup.total_out); + for (i = 0; i < zup.total_out; i++) { + fprintf (stderr, "%d : 0x%02hhx : 0x%02hhx\n", i, outbuf[i], srcbuf[i]); + assert (outbuf[i] == srcbuf[i]); + } + } + + exit (0); +} Index: trunk/Ohana/src/libfits/test/zlib.c =================================================================== --- trunk/Ohana/src/libfits/test/zlib.c (revision 38441) +++ trunk/Ohana/src/libfits/test/zlib.c (revision 38441) @@ -0,0 +1,210 @@ +# include +# include +# include +# include "tap_ohana.h" + +# define VERBOSE 0 +# define NSRC 10 +# define NTGT 20 +# define NOUT 125 + +# define SWAP_DBLE { \ + char tmp; \ + tmp = Pout[0]; Pout[0] = Pout[7]; Pout[7] = tmp; \ + tmp = Pout[1]; Pout[1] = Pout[6]; Pout[6] = tmp; \ + tmp = Pout[2]; Pout[2] = Pout[5]; Pout[5] = tmp; \ + tmp = Pout[3]; Pout[3] = Pout[4]; Pout[4] = tmp; } + +int main (int argc, char **argv) { + + plan_tests (7); + + diag ("libfits zlib.c tests"); + + diag ("zlib version: %s", ZLIB_VERSION); + diag ("zlib vernum: %x", ZLIB_VERNUM); + + int i, err; + double srcbuf[NSRC], tgtbuf[NTGT], outbuf[NOUT]; + + { + z_stream zdn; + zdn.next_in = (char *) srcbuf; + zdn.avail_in = NSRC * sizeof(double); + zdn.next_out = (char *) tgtbuf; + zdn.avail_out = NTGT * sizeof(double); + + zdn.zalloc = Z_NULL; + zdn.zfree = Z_NULL; + zdn.opaque = Z_NULL; + + for (i = 0; i < NSRC; i++) srcbuf[i] = 1.001*(i + 10); + memset (tgtbuf, 0, NTGT * sizeof(double)); + memset (outbuf, 0, NOUT * sizeof(double)); + + for (i = 0; i < NSRC; i++) { + char *Pout = (char *) &srcbuf[i]; + SWAP_DBLE; + } + + char *rawdata = (char *) srcbuf; + if (VERBOSE) fprintf (stderr, "inp: "); + for (i = 0; VERBOSE && (i < NSRC*8); i++) { + fprintf (stderr, "0x%02hhx ", rawdata[i]); + } + if (VERBOSE) fprintf (stderr, "\n"); + + // the '5' is the compression level: make this a user argument + err = deflateInit(&zdn, 5); + if (VERBOSE) fprintf (stderr, "inp buffers cmp 0: %d => %d => %d\n", zdn.avail_in, zdn.avail_out, (int) zdn.total_out); + ok (err == Z_OK, "result is Z_OK: %d vs %d", err, Z_OK); + + err = deflate(&zdn, Z_FINISH); + if (VERBOSE) fprintf (stderr, "out buffers cmp 1: %d => %d => %d\n", zdn.avail_in, zdn.avail_out, (int) zdn.total_out); + ok (err == Z_STREAM_END, "result is Z_STREAM_END: %d vs %d", err, Z_STREAM_END); + + // dump out the data which failed to uncompress (and the input data) + char *cmpdata = (char *) tgtbuf; + if (VERBOSE) fprintf (stderr, "out: "); + for (i = 0; VERBOSE && (i < zdn.total_out); i++) { + fprintf (stderr, "0x%02hhx ", cmpdata[i]); + } + if (VERBOSE) fprintf (stderr, "\n"); + + ok (zdn.total_out <= NTGT*sizeof(double), "output compressed size is good"); + + if (VERBOSE) fprintf (stderr, "result: %d bytes\n", (int) zdn.total_out); + for (i = 0; VERBOSE && (i < zdn.total_out / sizeof(double)); i++) { + fprintf (stderr, "%d : 0x%02hhx\n", i, tgtbuf[i]); + } + } + + { + z_stream zup; + zup.next_in = (char *) tgtbuf; + zup.avail_in = NTGT * sizeof(double); + zup.next_out = (char *) outbuf; + zup.avail_out = NOUT * sizeof(double); + + zup.zalloc = Z_NULL; + zup.zfree = Z_NULL; + zup.opaque = Z_NULL; + + // the '5' is the compression level: make this a user argument + err = inflateInit(&zup); + if (VERBOSE) fprintf (stderr, "out buffers unc 0: %d => %d => %d\n", zup.avail_in, zup.avail_out, (int) zup.total_out); + ok (err == Z_OK, "result is Z_OK: %d vs %d", err, Z_OK); + + err = inflate(&zup, Z_FINISH); + if (VERBOSE) fprintf (stderr, "out buffers unc 1: %d => %d => %d\n", zup.avail_in, zup.avail_out, (int) zup.total_out); + ok (err == Z_STREAM_END, "result is Z_STREAM_END: %d vs %d", err, Z_STREAM_END); + + char *outdata = (char *) outbuf; + if (VERBOSE) fprintf (stderr, "unc: "); + for (i = 0; VERBOSE && (i < zup.total_out); i++) { + fprintf (stderr, "0x%02hhx ", outdata[i]); + } + if (VERBOSE) fprintf (stderr, "\n"); + + ok (zup.total_out <= NTGT*sizeof(double), "uncompressed size is good"); + + int Nbad = 0; + for (i = 0; 0 && i < zup.total_out / sizeof(double); i++) { + if (outbuf[i] != srcbuf[i]) Nbad++; + } + ok (Nbad == 0, "no mismatched bytes"); + } + + /**** TEST 2 : use gzip headers ****/ + { + z_stream zdn; + zdn.next_in = (char *) srcbuf; + zdn.avail_in = NSRC * sizeof(double); + zdn.next_out = (char *) tgtbuf; + zdn.avail_out = NTGT * sizeof(double); + + zdn.zalloc = Z_NULL; + zdn.zfree = Z_NULL; + zdn.opaque = Z_NULL; + + for (i = 0; i < NSRC; i++) srcbuf[i] = 1.001*(i + 10); + memset (tgtbuf, 0, NTGT * sizeof(double)); + memset (outbuf, 0, NOUT * sizeof(double)); + + for (i = 0; i < NSRC; i++) { + char *Pout = (char *) &srcbuf[i]; + SWAP_DBLE; + } + + char *rawdata = (char *) srcbuf; + if (VERBOSE) fprintf (stderr, "inp: "); + for (i = 0; VERBOSE && (i < NSRC*8); i++) { + fprintf (stderr, "0x%02hhx ", rawdata[i]); + } + if (VERBOSE) fprintf (stderr, "\n"); + + // the '5' is the compression level: make this a user argument + // windowBits = 31 = (15 + 16) = (2^15 window bits) + (create a gzip stream) + err = deflateInit2(&zdn, 1, Z_DEFLATED, 31, 8, Z_DEFAULT_STRATEGY); + if (VERBOSE) fprintf (stderr, "inp buffers cmp 0: %d => %d => %d\n", zdn.avail_in, zdn.avail_out, (int) zdn.total_out); + ok (err == Z_OK, "result is Z_OK: %d vs %d", err, Z_OK); + + err = deflate(&zdn, Z_FINISH); + if (VERBOSE) fprintf (stderr, "out buffers cmp 1: %d => %d => %d\n", zdn.avail_in, zdn.avail_out, (int) zdn.total_out); + ok (err == Z_STREAM_END, "result is Z_STREAM_END: %d vs %d", err, Z_STREAM_END); + + // dump out the data which failed to uncompress (and the input data) + char *cmpdata = (char *) tgtbuf; + if (VERBOSE) fprintf (stderr, "out: "); + for (i = 0; VERBOSE && (i < zdn.total_out); i++) { + fprintf (stderr, "0x%02hhx ", cmpdata[i]); + } + if (VERBOSE) fprintf (stderr, "\n"); + + ok (zdn.total_out <= NTGT*sizeof(double), "output compressed size is good"); + + if (VERBOSE) fprintf (stderr, "result: %d bytes\n", (int) zdn.total_out); + for (i = 0; VERBOSE && (i < zdn.total_out / sizeof(double)); i++) { + fprintf (stderr, "%d : 0x%02hhx\n", i, tgtbuf[i]); + } + } + + { + z_stream zup; + zup.next_in = (char *) tgtbuf; + zup.avail_in = NTGT * sizeof(double); + zup.next_out = (char *) outbuf; + zup.avail_out = NOUT * sizeof(double); + + zup.zalloc = Z_NULL; + zup.zfree = Z_NULL; + zup.opaque = Z_NULL; + + // a value of 32 tells inflate to use the window size used in the compression AND to look for both zlib and gzip headers + // windowBits = 47 = (15 + 32) = (2^15 window bits) + (test for either gzip or zlib streams) + err = inflateInit2(&zup, 47); + if (VERBOSE) fprintf (stderr, "out buffers unc 0: %d => %d => %d\n", zup.avail_in, zup.avail_out, (int) zup.total_out); + ok (err == Z_OK, "inflateInit2 result is Z_OK: %d vs %d", err, Z_OK); + + err = inflate(&zup, Z_FINISH); + if (VERBOSE) fprintf (stderr, "out buffers unc 1: %d => %d => %d\n", zup.avail_in, zup.avail_out, (int) zup.total_out); + ok (err == Z_STREAM_END, "inflate result is Z_STREAM_END: %d vs %d", err, Z_STREAM_END); + + char *outdata = (char *) outbuf; + if (VERBOSE) fprintf (stderr, "unc: "); + for (i = 0; VERBOSE && (i < zup.total_out); i++) { + fprintf (stderr, "0x%02hhx ", outdata[i]); + } + if (VERBOSE) fprintf (stderr, "\n"); + + ok (zup.total_out <= NTGT*sizeof(double), "uncompressed size is good"); + + int Nbad = 0; + for (i = 0; 0 && i < zup.total_out / sizeof(double); i++) { + if (outbuf[i] != srcbuf[i]) Nbad++; + } + ok (Nbad == 0, "no mismatched bytes"); + } + + return exit_status(); +}