Index: trunk/Ohana/src/libfits/Makefile =================================================================== --- trunk/Ohana/src/libfits/Makefile (revision 15240) +++ trunk/Ohana/src/libfits/Makefile (revision 15487) @@ -10,4 +10,5 @@ MATR = $(HOME)/matrix TABL = $(HOME)/table +EXT = $(HOME)/extern INC = $(HOME)/include MAN = $(HOME)/doc @@ -15,6 +16,6 @@ # programs may add their own internal requirements here -FULL_CFLAGS = $(BASE_CFLAGS) -fPIC -FULL_CPPFLAGS = $(BASE_CPPFLAGS) +FULL_CFLAGS = $(BASE_CFLAGS) -fPIC -Wall -Werror +FULL_CPPFLAGS = $(BASE_CPPFLAGS) -I$(EXT) FULL_LDFLAGS = $(BASE_LDFLAGS) @@ -32,5 +33,5 @@ $(HEAD)/F_copy_H.$(ARCH).o $(HEAD)/F_delete.$(ARCH).o \ $(HEAD)/F_read_XH.$(ARCH).o $(HEAD)/F_init_H.$(ARCH).o \ -$(HEAD)/version.$(ARCH).o +$(HEAD)/F_convert_H.$(ARCH).o $(HEAD)/version.$(ARCH).o MATRIX_OBJ = \ @@ -43,5 +44,6 @@ $(MATR)/F_convert_format.$(ARCH).o $(MATR)/F_read_segment.$(ARCH).o \ $(MATR)/F_read_portion.$(ARCH).o $(MATR)/F_load_M.$(ARCH).o \ -$(MATR)/F_matrix.$(ARCH).o +$(MATR)/F_matrix.$(ARCH).o $(MATR)/F_compress_M.$(ARCH).o \ +$(MATR)/F_uncompress_data.$(ARCH).o TABLE_OBJ = \ @@ -51,7 +53,16 @@ $(TABL)/F_define_column.$(ARCH).o $(TABL)/F_table_format.$(ARCH).o \ $(TABL)/F_set_column.$(ARCH).o $(TABL)/F_get_column.$(ARCH).o \ -$(TABL)/F_table_row.$(ARCH).o $(TABL)/F_free_T.$(ARCH).o +$(TABL)/F_table_row.$(ARCH).o $(TABL)/F_free_T.$(ARCH).o \ +$(TABL)/F_table_varlength.$(ARCH).o -OBJS = $(HEADER_OBJ) $(MATRIX_OBJ) $(TABLE_OBJ) +EXTERN_OBJ = \ +$(EXT)/fits_hcompress.$(ARCH).o \ +$(EXT)/fits_hdecompress.$(ARCH).o \ +$(EXT)/pliocomp.$(ARCH).o \ +$(EXT)/ricecomp.$(ARCH).o + +OBJS = $(HEADER_OBJ) $(MATRIX_OBJ) $(TABLE_OBJ) $(EXTERN_OBJ) + +extern: $(EXTERN_OBJ) $(OBJS): $(INCS) Index: trunk/Ohana/src/libfits/doc/notes-compress.txt =================================================================== --- trunk/Ohana/src/libfits/doc/notes-compress.txt (revision 15240) +++ trunk/Ohana/src/libfits/doc/notes-compress.txt (revision 15487) @@ -1,2 +1,16 @@ + +TFORMn : rPt(e_max) +r is 0, 1, or absent +P is literal +t is datatype character for binary tables +e_max is the largest record in the column + +data column actually contains: 2 32bit values: +Length + Offset + +THEAP : offset from start of data to start of heap +(note that the heap starts at NAXIS1 x NAXIS2, not Nx2880 bytes). + +total size of heap segment (gap + heap data) is PCOUNT I need to implement (at least) support for read of compressed images. @@ -16,5 +30,5 @@ level ones as needed. -* possible needed APIs: +* possible needed APIs & upgrades: - gfits_image_is_compressed (Header *header) @@ -30,3 +44,11 @@ compress the given image returning the table data + o include PCOUNT in data area size (gfits_matrix_size) + o read heap when reading table (gfits_fread_table_data) + + - void *gfits_varlength_column_pointer (FTable *ftable, VarLengthColumn *column, int row); + + o int gfits_varlength_column_define (FTable *ftable, VarLengthColumn *def, int column); + + Index: trunk/Ohana/src/libfits/extern/fits_hcompress.c =================================================================== --- trunk/Ohana/src/libfits/extern/fits_hcompress.c (revision 15487) +++ trunk/Ohana/src/libfits/extern/fits_hcompress.c (revision 15487) @@ -0,0 +1,1688 @@ +/* ######################################################################### +These routines to apply the H-compress compression algorithm to a 2-D Fits +image were written by R. White at the STScI and were obtained from the STScI at +http://www.stsci.edu/software/hcompress.html + +This source file is a concatination of the following sources files in the +original distribution + htrans.c + digitize.c + encode.c + qwrite.c + doencode.c + bit_output.c + qtree_encode.c + +The following modifications have been made to the original code: + + - commented out redundant "include" statements + - added the noutchar global variable + - changed all the 'extern' declarations to 'static', since all the routines are in + the same source file + - changed the first parameter in encode (and in lower level routines from a file stream + to a char array + - modifid the encode routine to return the size of the compressed array of bytes + - changed calls to printf and perror to call the CFITSIO ffpmsg routine + - modified the mywrite routine, and lower level byte writing routines, to copy + the output bytes to a char array, instead of writing them to a file stream + - replace "exit" statements with "return" statements + - changed the function declarations to the more modern ANSI C style + + ############################################################################ */ + +#include +#include +#include +#include + +// EAM : dropping fitsio prototypes +# define LONGLONG long long +# define DATA_COMPRESSION_ERR 413 /* error in imcompress routines */ +# define DATA_DECOMPRESSION_ERR 414 /* error in imcompress routines */ + +static long noutchar; +static long noutmax; + +static int htrans(int a[],int nx,int ny); +static void digitize(int a[], int nx, int ny, int scale); +static int encode(char *outfile, long *nlen, int a[], int nx, int ny, int scale); +static void shuffle(int a[], int n, int n2, int tmp[]); + +static int htrans64(LONGLONG a[],int nx,int ny); +static void digitize64(LONGLONG a[], int nx, int ny, int scale); +static int encode64(char *outfile, long *nlen, LONGLONG a[], int nx, int ny, int scale); +static void shuffle64(LONGLONG a[], int n, int n2, LONGLONG tmp[]); + +static void writeint(char *outfile, int a); +static void writelonglong(char *outfile, LONGLONG a); +static int qwrite(char *outfile, char *a, int n); +static int doencode(char *outfile, int a[], int nx, int ny, unsigned char nbitplanes[3]); +static int doencode64(char *outfile, LONGLONG a[], int nx, int ny, unsigned char nbitplanes[3]); +static int mywrite(char *file, char buffer[], int n); + +static int qtree_encode(char *outfile, int a[], int n, int nqx, int nqy, int nbitplanes); +static int qtree_encode64(char *outfile, LONGLONG a[], int n, int nqx, int nqy, int nbitplanes); +static void start_outputing_bits(); +static void done_outputing_bits(char *outfile); +static void output_nbits(char *outfile, int bits, int n); + +static void qtree_onebit(int a[], int n, int nx, int ny, unsigned char b[], int bit); +static void qtree_onebit64(LONGLONG a[], int n, int nx, int ny, unsigned char b[], int bit); +static void qtree_reduce(unsigned char a[], int n, int nx, int ny, unsigned char b[]); +static int bufcopy(unsigned char a[], int n, unsigned char buffer[], int *b, int bmax); +static void write_bdirect(char *outfile, int a[], int n,int nqx, int nqy, unsigned char scratch[], int bit); +static void write_bdirect64(char *outfile, LONGLONG a[], int n,int nqx, int nqy, unsigned char scratch[], int bit); + +#define output_nybble(outfile,c) output_nbits(outfile,c,4) +#define output_huffman(outfile,c) output_nbits(outfile,code[c],ncode[c]) + +/* ---------------------------------------------------------------------- */ +int fits_hcompress(int *a, int ny, int nx, int scale, char *output, + long *nbytes, int *status) +{ + /* + compress the input image using the H-compress algorithm + + a - input image array + nx - size of X axis of image + ny - size of Y axis of image + scale - quantization scale factor. Larger values results in more (lossy) compression + scale = 0 does lossless compression + output - pre-allocated array to hold the output compressed stream of bytes + nbyts - input value = size of the output buffer; + returned value = size of the compressed byte stream, in bytes + + NOTE: the nx and ny dimensions as defined within this code are reversed from + the usual FITS notation. ny is the fastest varying dimension, which is + usually considered the X axis in the FITS image display + + */ + + int stat; + + if (*status > 0) return(*status); + + /* H-transform */ + stat = htrans(a, nx, ny); + if (stat) { + *status = stat; + return(*status); + } + + /* digitize */ + digitize(a, nx, ny, scale); + + /* encode and write to output array */ + + noutmax = *nbytes; /* input value is the allocated size of the array */ + *nbytes = 0; /* reset */ + + stat = encode(output, nbytes, a, nx, ny, scale); + + *status = stat; + return(*status); +} +/* ---------------------------------------------------------------------- */ +int fits_hcompress64(LONGLONG *a, int ny, int nx, int scale, char *output, + long *nbytes, int *status) +{ + /* + compress the input image using the H-compress algorithm + + a - input image array + nx - size of X axis of image + ny - size of Y axis of image + scale - quantization scale factor. Larger values results in more (lossy) compression + scale = 0 does lossless compression + output - pre-allocated array to hold the output compressed stream of bytes + nbyts - size of the compressed byte stream, in bytes + + NOTE: the nx and ny dimensions as defined within this code are reversed from + the usual FITS notation. ny is the fastest varying dimension, which is + usually considered the X axis in the FITS image display + + */ + + int stat; + + if (*status > 0) return(*status); + + /* H-transform */ + stat = htrans64(a, nx, ny); + if (stat) { + *status = stat; + return(*status); + } + + /* digitize */ + digitize64(a, nx, ny, scale); + + /* encode and write to output array */ + noutmax = *nbytes; /* input value is the allocated size of the array */ + *nbytes = 0; /* reset */ + + stat = encode64(output, nbytes, a, nx, ny, scale); + + *status = stat; + return(*status); +} + + +/* Copyright (c) 1993 Association of Universities for Research + * in Astronomy. All rights reserved. Produced under National + * Aeronautics and Space Administration Contract No. NAS5-26555. + */ +/* htrans.c H-transform of NX x NY integer image + * + * Programmer: R. White Date: 11 May 1992 + */ + +/* ######################################################################### */ +static int htrans(int a[],int nx,int ny) +{ +int nmax, log2n, h0, hx, hy, hc, nxtop, nytop, i, j, k; +int oddx, oddy; +int shift, mask, mask2, prnd, prnd2, nrnd2; +int s10, s00; +int *tmp; + + /* + * log2n is log2 of max(nx,ny) rounded up to next power of 2 + */ + nmax = (nx>ny) ? nx : ny; + log2n = (int) (log((float) nmax)/log(2.0)+0.5); + if ( nmax > (1<> shift; + hx = (a[s10+1] + a[s10] - a[s00+1] - a[s00]) >> shift; + hy = (a[s10+1] - a[s10] + a[s00+1] - a[s00]) >> shift; + hc = (a[s10+1] - a[s10] - a[s00+1] + a[s00]) >> shift; + /* + * Throw away the 2 bottom bits of h0, bottom bit of hx,hy. + * To get rounding to be same for positive and negative + * numbers, nrnd2 = prnd2 - 1. + */ + a[s10+1] = hc; + a[s10 ] = ( (hx>=0) ? (hx+prnd) : hx ) & mask ; + a[s00+1] = ( (hy>=0) ? (hy+prnd) : hy ) & mask ; + a[s00 ] = ( (h0>=0) ? (h0+prnd2) : (h0+nrnd2) ) & mask2; + s00 += 2; + s10 += 2; + } + if (oddy) { + /* + * do last element in row if row length is odd + * s00+1, s10+1 are off edge + */ + h0 = (a[s10] + a[s00]) << (1-shift); + hx = (a[s10] - a[s00]) << (1-shift); + a[s10 ] = ( (hx>=0) ? (hx+prnd) : hx ) & mask ; + a[s00 ] = ( (h0>=0) ? (h0+prnd2) : (h0+nrnd2) ) & mask2; + s00 += 1; + s10 += 1; + } + } + if (oddx) { + /* + * do last row if column length is odd + * s10, s10+1 are off edge + */ + s00 = i*ny; + for (j = 0; j=0) ? (hy+prnd) : hy ) & mask ; + a[s00 ] = ( (h0>=0) ? (h0+prnd2) : (h0+nrnd2) ) & mask2; + s00 += 2; + } + if (oddy) { + /* + * do corner element if both row and column lengths are odd + * s00+1, s10, s10+1 are off edge + */ + h0 = a[s00] << (2-shift); + a[s00 ] = ( (h0>=0) ? (h0+prnd2) : (h0+nrnd2) ) & mask2; + } + } + /* + * now shuffle in each dimension to group coefficients by order + */ + for (i = 0; i>1; + nytop = (nytop+1)>>1; + /* + * divisor doubles after first reduction + */ + shift = 1; + /* + * masks, rounding values double after each iteration + */ + mask = mask2; + prnd = prnd2; + mask2 = mask2 << 1; + prnd2 = prnd2 << 1; + nrnd2 = prnd2 - 1; + } + free(tmp); + return(0); +} +/* ######################################################################### */ + +static int htrans64(LONGLONG a[],int nx,int ny) +{ +int nmax, log2n, nxtop, nytop, i, j, k; +int oddx, oddy; +int shift; +int s10, s00; +LONGLONG h0, hx, hy, hc, prnd, prnd2, nrnd2, mask, mask2; +LONGLONG *tmp; + + /* + * log2n is log2 of max(nx,ny) rounded up to next power of 2 + */ + nmax = (nx>ny) ? nx : ny; + log2n = (int) (log((float) nmax)/log(2.0)+0.5); + if ( nmax > (1<> shift; + hx = (a[s10+1] + a[s10] - a[s00+1] - a[s00]) >> shift; + hy = (a[s10+1] - a[s10] + a[s00+1] - a[s00]) >> shift; + hc = (a[s10+1] - a[s10] - a[s00+1] + a[s00]) >> shift; + /* + * Throw away the 2 bottom bits of h0, bottom bit of hx,hy. + * To get rounding to be same for positive and negative + * numbers, nrnd2 = prnd2 - 1. + */ + a[s10+1] = hc; + a[s10 ] = ( (hx>=0) ? (hx+prnd) : hx ) & mask ; + a[s00+1] = ( (hy>=0) ? (hy+prnd) : hy ) & mask ; + a[s00 ] = ( (h0>=0) ? (h0+prnd2) : (h0+nrnd2) ) & mask2; + s00 += 2; + s10 += 2; + } + if (oddy) { + /* + * do last element in row if row length is odd + * s00+1, s10+1 are off edge + */ + h0 = (a[s10] + a[s00]) << (1-shift); + hx = (a[s10] - a[s00]) << (1-shift); + a[s10 ] = ( (hx>=0) ? (hx+prnd) : hx ) & mask ; + a[s00 ] = ( (h0>=0) ? (h0+prnd2) : (h0+nrnd2) ) & mask2; + s00 += 1; + s10 += 1; + } + } + if (oddx) { + /* + * do last row if column length is odd + * s10, s10+1 are off edge + */ + s00 = i*ny; + for (j = 0; j=0) ? (hy+prnd) : hy ) & mask ; + a[s00 ] = ( (h0>=0) ? (h0+prnd2) : (h0+nrnd2) ) & mask2; + s00 += 2; + } + if (oddy) { + /* + * do corner element if both row and column lengths are odd + * s00+1, s10, s10+1 are off edge + */ + h0 = a[s00] << (2-shift); + a[s00 ] = ( (h0>=0) ? (h0+prnd2) : (h0+nrnd2) ) & mask2; + } + } + /* + * now shuffle in each dimension to group coefficients by order + */ + for (i = 0; i>1; + nytop = (nytop+1)>>1; + /* + * divisor doubles after first reduction + */ + shift = 1; + /* + * masks, rounding values double after each iteration + */ + mask = mask2; + prnd = prnd2; + mask2 = mask2 << 1; + prnd2 = prnd2 << 1; + nrnd2 = prnd2 - 1; + } + free(tmp); + return(0); +} + +/* ######################################################################### */ +static void +shuffle(int a[], int n, int n2, int tmp[]) +{ + +/* +int a[]; array to shuffle +int n; number of elements to shuffle +int n2; second dimension +int tmp[]; scratch storage +*/ + +int i; +int *p1, *p2, *pt; + + /* + * copy odd elements to tmp + */ + pt = tmp; + p1 = &a[n2]; + for (i=1; i < n; i += 2) { + *pt = *p1; + pt += 1; + p1 += (n2+n2); + } + /* + * compress even elements into first half of A + */ + p1 = &a[n2]; + p2 = &a[n2+n2]; + for (i=2; i0) ? (*p+d) : (*p-d))/scale; +} + +/* ######################################################################### */ +static void +digitize64(LONGLONG a[], int nx, int ny, int scale) +{ +LONGLONG d, *p, scale64; + + /* + * round to multiple of scale + */ + if (scale <= 1) return; + d=(scale+1)/2-1; + scale64 = scale; /* use a 64-bit int for efficiency in the big loop */ + + for (p=a; p <= &a[nx*ny-1]; p++) *p = ((*p>0) ? (*p+d) : (*p-d))/scale64; +} +/* ######################################################################### */ +/* ######################################################################### */ +/* Copyright (c) 1993 Association of Universities for Research + * in Astronomy. All rights reserved. Produced under National + * Aeronautics and Space Administration Contract No. NAS5-26555. + */ +/* encode.c encode H-transform and write to outfile + * + * Programmer: R. White Date: 2 February 1994 + */ + +static char code_magic[2] = { (char)0xDD, (char)0x99 }; + + +/* ######################################################################### */ +static int encode(char *outfile, long *nlength, int a[], int nx, int ny, int scale) +{ + +/* FILE *outfile; - change outfile to a char array */ +/* + long * nlength returned length (in bytes) of the encoded array) + int a[]; input H-transform array (nx,ny) + int nx,ny; size of H-transform array + int scale; scale factor for digitization +*/ +int nel, nx2, ny2, i, j, k, q, vmax[3], nsign, bits_to_go; +unsigned char nbitplanes[3]; +unsigned char *signbits; +int stat = 0; + + noutchar = 0; /* initialize the number of compressed bytes that have been written */ + nel = nx*ny; + /* + * write magic value + */ + qwrite(outfile, code_magic, sizeof(code_magic)); + writeint(outfile, nx); /* size of image */ + writeint(outfile, ny); + writeint(outfile, scale); /* scale factor for digitization */ + /* + * write first value of A (sum of all pixels -- the only value + * which does not compress well) + */ + writelonglong(outfile, (LONGLONG) a[0]); + + a[0] = 0; + /* + * allocate array for sign bits and save values, 8 per byte + */ + signbits = (unsigned char *) malloc((nel+7)/8); + if (signbits == (unsigned char *) NULL) { + // EAM : dropping Pence CFITSIO dependencies + // XXX : add an error message function to gfits + fprintf(stderr, "encode: insufficient memory\n"); + return(DATA_COMPRESSION_ERR); + } + nsign = 0; + bits_to_go = 8; + signbits[0] = 0; + for (i=0; i 0) { + /* + * positive element, put zero at end of buffer + */ + signbits[nsign] <<= 1; + bits_to_go -= 1; + } else if (a[i] < 0) { + /* + * negative element, shift in a one + */ + signbits[nsign] <<= 1; + signbits[nsign] |= 1; + bits_to_go -= 1; + /* + * replace a by absolute value + */ + a[i] = -a[i]; + } + if (bits_to_go == 0) { + /* + * filled up this byte, go to the next one + */ + bits_to_go = 8; + nsign += 1; + signbits[nsign] = 0; + } + } + if (bits_to_go != 8) { + /* + * some bits in last element + * move bits in last byte to bottom and increment nsign + */ + signbits[nsign] <<= bits_to_go; + nsign += 1; + } + + /* + * calculate number of bit planes for 3 quadrants + * + * quadrant 0=bottom left, 1=bottom right or top left, 2=top right, + */ + for (q=0; q<3; q++) { + vmax[q] = 0; + } + /* + * get maximum absolute value in each quadrant + */ + nx2 = (nx+1)/2; + ny2 = (ny+1)/2; + j=0; /* column counter */ + k=0; /* row counter */ + for (i=0; i=ny2) + (k>=nx2); + if (vmax[q] < a[i]) vmax[q] = a[i]; + if (++j >= ny) { + j = 0; + k += 1; + } + } + /* + * now calculate number of bits for each quadrant + */ + for (q = 0; q < 3; q++) { + nbitplanes[q] = (int) (log((float) (vmax[q]+1))/log(2.0)+0.5); + if ( (vmax[q]+1) > (1< 0) { + + if ( 0 == qwrite(outfile, (char *) signbits, nsign)) { + free(signbits); + *nlength = noutchar; + fprintf (stderr, "encode: output buffer too small"); + return(DATA_COMPRESSION_ERR); + } + } + + free(signbits); + *nlength = noutchar; + + if (noutchar >= noutmax) { + fprintf (stderr, "encode: output buffer too small"); + return(DATA_COMPRESSION_ERR); + } + + return(stat); +} +/* ######################################################################### */ +static int encode64(char *outfile, long *nlength, LONGLONG a[], int nx, int ny, int scale) +{ + +/* FILE *outfile; - change outfile to a char array */ +/* + long * nlength returned length (in bytes) of the encoded array) + LONGLONG a[]; input H-transform array (nx,ny) + int nx,ny; size of H-transform array + int scale; scale factor for digitization +*/ +int nel, nx2, ny2, i, j, k, q, nsign, bits_to_go; +LONGLONG vmax[3]; +unsigned char nbitplanes[3]; +unsigned char *signbits; +int stat = 0; + + noutchar = 0; /* initialize the number of compressed bytes that have been written */ + nel = nx*ny; + /* + * write magic value + */ + qwrite(outfile, code_magic, sizeof(code_magic)); + writeint(outfile, nx); /* size of image */ + writeint(outfile, ny); + writeint(outfile, scale); /* scale factor for digitization */ + /* + * write first value of A (sum of all pixels -- the only value + * which does not compress well) + */ + + writelonglong(outfile, a[0]); + a[0] = 0; + /* + * allocate array for sign bits and save values, 8 per byte + */ + signbits = (unsigned char *) malloc((nel+7)/8); + if (signbits == (unsigned char *) NULL) { + fprintf (stderr, "encode64: insufficient memory"); + return(DATA_COMPRESSION_ERR); + } + nsign = 0; + bits_to_go = 8; + signbits[0] = 0; + for (i=0; i 0) { + /* + * positive element, put zero at end of buffer + */ + signbits[nsign] <<= 1; + bits_to_go -= 1; + } else if (a[i] < 0) { + /* + * negative element, shift in a one + */ + signbits[nsign] <<= 1; + signbits[nsign] |= 1; + bits_to_go -= 1; + /* + * replace a by absolute value + */ + a[i] = -a[i]; + } + if (bits_to_go == 0) { + /* + * filled up this byte, go to the next one + */ + bits_to_go = 8; + nsign += 1; + signbits[nsign] = 0; + } + } + if (bits_to_go != 8) { + /* + * some bits in last element + * move bits in last byte to bottom and increment nsign + */ + signbits[nsign] <<= bits_to_go; + nsign += 1; + } + /* + * calculate number of bit planes for 3 quadrants + * + * quadrant 0=bottom left, 1=bottom right or top left, 2=top right, + */ + for (q=0; q<3; q++) { + vmax[q] = 0; + } + /* + * get maximum absolute value in each quadrant + */ + nx2 = (nx+1)/2; + ny2 = (ny+1)/2; + + i = 0; + for (k=0; k=ny2) + (k>=nx2); + vmax[q] |= a[i]; + i++; + } + } + + /* + * now calculate number of bits for each quadrant + */ + + /* this is a more efficient way to do this, */ + + + for (q = 0; q < 3; q++) { + for (nbitplanes[q] = 0; vmax[q]>0; vmax[q] = vmax[q]>>1, nbitplanes[q]++) ; + } + + +/* + for (q = 0; q < 3; q++) { + nbitplanes[q] = log((float) (vmax[q]+1))/log(2.0)+0.5; + if ( (vmax[q]+1) > (((LONGLONG) 1)< 0) { + + if ( 0 == qwrite(outfile, (char *) signbits, nsign)) { + free(signbits); + *nlength = noutchar; + fprintf (stderr, "encode: output buffer too small"); + return(DATA_COMPRESSION_ERR); + } + } + + free(signbits); + + *nlength = noutchar; + + if (noutchar >= noutmax) { + fprintf (stderr, "encode64: output buffer too small"); + return(DATA_COMPRESSION_ERR); + } + + return(stat); +} +/* ######################################################################### */ +/* ######################################################################### */ +/* Copyright (c) 1993 Association of Universities for Research + * in Astronomy. All rights reserved. Produced under National + * Aeronautics and Space Administration Contract No. NAS5-26555. + */ +/* qwrite.c Write binary data + * + * Programmer: R. White Date: 11 March 1991 + */ + +/* ######################################################################### */ +static void +writeint(char *outfile, int a) +{ +int i; +unsigned char b[4]; + + /* Write integer A one byte at a time to outfile. + * + * This is portable from Vax to Sun since it eliminates the + * need for byte-swapping. + */ + for (i=3; i>=0; i--) { + b[i] = a & 0x000000ff; + a >>= 8; + } + for (i=0; i<4; i++) qwrite(outfile, (char *) &b[i],1); +} + +/* ######################################################################### */ +static void +writelonglong(char *outfile, LONGLONG a) +{ +int i; +unsigned char b[8]; + + /* Write integer A one byte at a time to outfile. + * + * This is portable from Vax to Sun since it eliminates the + * need for byte-swapping. + */ + for (i=7; i>=0; i--) { + b[i] = (unsigned char) (a & 0x000000ff); + a >>= 8; + } + for (i=0; i<8; i++) qwrite(outfile, (char *) &b[i],1); +} + +/* ######################################################################### */ +static int qwrite(char *outfile, char *a, int n) +{ +int nwrite; + + nwrite = mywrite(outfile, a, n); + return(nwrite); /* added by WDP to prevent compiler warning */ +} + +/* ######################################################################### */ +static int +mywrite(char *file, char buffer[], int n) +{ + /* + * write n bytes from buffer into file + * returns number of bytes read (=n) if successful, <=0 if not + */ + + if (noutchar + n > noutmax) return(0); /* buffer overflow */ + + memcpy(&file[noutchar], buffer, n); + noutchar += n; + + return(n); +} +/* ######################################################################### */ +/* ######################################################################### */ +/* Copyright (c) 1993 Association of Universities for Research + * in Astronomy. All rights reserved. Produced under National + * Aeronautics and Space Administration Contract No. NAS5-26555. + */ +/* doencode.c Encode 2-D array and write stream of characters on outfile + * + * This version assumes that A is positive. + * + * Programmer: R. White Date: 7 May 1991 + */ + +/* ######################################################################### */ +static int +doencode(char *outfile, int a[], int nx, int ny, unsigned char nbitplanes[3]) +{ +/* char *outfile; output data stream +int a[]; Array of values to encode +int nx,ny; Array dimensions [nx][ny] +unsigned char nbitplanes[3]; Number of bit planes in quadrants +*/ + +int nx2, ny2, stat = 0; + + nx2 = (nx+1)/2; + ny2 = (ny+1)/2; + /* + * Initialize bit output + */ + start_outputing_bits(); + /* + * write out the bit planes for each quadrant + */ + stat = qtree_encode(outfile, &a[0], ny, nx2, ny2, nbitplanes[0]); + + if (!stat) + stat = qtree_encode(outfile, &a[ny2], ny, nx2, ny/2, nbitplanes[1]); + + if (!stat) + stat = qtree_encode(outfile, &a[ny*nx2], ny, nx/2, ny2, nbitplanes[1]); + + if (!stat) + stat = qtree_encode(outfile, &a[ny*nx2+ny2], ny, nx/2, ny/2, nbitplanes[2]); + /* + * Add zero as an EOF symbol + */ + output_nybble(outfile, 0); + done_outputing_bits(outfile); + + return(stat); +} +/* ######################################################################### */ +static int +doencode64(char *outfile, LONGLONG a[], int nx, int ny, unsigned char nbitplanes[3]) +{ +/* char *outfile; output data stream +LONGLONG a[]; Array of values to encode +int nx,ny; Array dimensions [nx][ny] +unsigned char nbitplanes[3]; Number of bit planes in quadrants +*/ + +int nx2, ny2, stat = 0; + + nx2 = (nx+1)/2; + ny2 = (ny+1)/2; + /* + * Initialize bit output + */ + start_outputing_bits(); + /* + * write out the bit planes for each quadrant + */ + stat = qtree_encode64(outfile, &a[0], ny, nx2, ny2, nbitplanes[0]); + + if (!stat) + stat = qtree_encode64(outfile, &a[ny2], ny, nx2, ny/2, nbitplanes[1]); + + if (!stat) + stat = qtree_encode64(outfile, &a[ny*nx2], ny, nx/2, ny2, nbitplanes[1]); + + if (!stat) + stat = qtree_encode64(outfile, &a[ny*nx2+ny2], ny, nx/2, ny/2, nbitplanes[2]); + /* + * Add zero as an EOF symbol + */ + output_nybble(outfile, 0); + done_outputing_bits(outfile); + + return(stat); +} +/* ######################################################################### */ +/* ######################################################################### */ +/* Copyright (c) 1993 Association of Universities for Research + * in Astronomy. All rights reserved. Produced under National + * Aeronautics and Space Administration Contract No. NAS5-26555. + */ +/* BIT OUTPUT ROUTINES */ + + +static LONGLONG bitcount; + +/* THE BIT BUFFER */ + +static int buffer2; /* Bits buffered for output */ +static int bits_to_go2; /* Number of bits free in buffer */ + + +/* ######################################################################### */ +/* INITIALIZE FOR BIT OUTPUT */ + +static void +start_outputing_bits() +{ + buffer2 = 0; /* Buffer is empty to start */ + bits_to_go2 = 8; /* with */ + bitcount = 0; +} + +/* ######################################################################### */ +/* OUTPUT N BITS (N must be <= 8) */ + +static void +output_nbits(char *outfile, int bits, int n) +{ + /* + * insert bits at end of buffer + */ + buffer2 <<= n; + buffer2 |= ( bits & ((1<>(-bits_to_go2)) & 0xff,outfile); */ + + outfile[noutchar] = ((buffer2>>(-bits_to_go2)) & 0xff); + + if (noutchar < noutmax) noutchar++; + + bits_to_go2 += 8; + } + bitcount += n; +} + + +/* ######################################################################### */ +/* FLUSH OUT THE LAST BITS */ + +static void +done_outputing_bits(char *outfile) +{ + if(bits_to_go2 < 8) { +/* putc(buffer2<nqy) ? nqx : nqy; + log2n = (int) (log((float) nqmax)/log(2.0)+0.5); + if (nqmax > (1<= 0; bit--) { + /* + * initial bit buffer + */ + b = 0; + bitbuffer = 0; + bits_to_go3 = 0; + /* + * on first pass copy A to scratch array + */ + qtree_onebit(a,n,nqx,nqy,scratch,bit); + nx = (nqx+1)>>1; + ny = (nqy+1)>>1; + /* + * copy non-zero values to output buffer, which will be written + * in reverse order + */ + if (bufcopy(scratch,nx*ny,buffer,&b,bmax)) { + /* + * quadtree is expanding data, + * change warning code and just fill buffer with bit-map + */ + write_bdirect(outfile,a,n,nqx,nqy,scratch,bit); + goto bitplane_done; + } + /* + * do log2n reductions + */ + for (k = 1; k>1; + ny = (ny+1)>>1; + if (bufcopy(scratch,nx*ny,buffer,&b,bmax)) { + write_bdirect(outfile,a,n,nqx,nqy,scratch,bit); + goto bitplane_done; + } + } + /* + * OK, we've got the code in buffer + * Write quadtree warning code, then write buffer in reverse order + */ + output_nybble(outfile,0xF); + if (b==0) { + if (bits_to_go3>0) { + /* + * put out the last few bits + */ + output_nbits(outfile, bitbuffer & ((1<0) { + /* + * put out the last few bits + */ + output_nbits(outfile, bitbuffer & ((1<=0; i--) { + output_nbits(outfile,buffer[i],8); + } + } + bitplane_done: ; + } + free(buffer); + free(scratch); + return(0); +} +/* ######################################################################### */ +static int +qtree_encode64(char *outfile, LONGLONG a[], int n, int nqx, int nqy, int nbitplanes) +{ + +/* +LONGLONG a[]; +int n; physical dimension of row in a +int nqx; length of row +int nqy; length of column (<=n) +int nbitplanes; number of bit planes to output +*/ + +int log2n, i, k, bit, b, nqmax, nqx2, nqy2, nx, ny; +int bmax; /* this potentially needs to be made a 64-bit int to support large arrays */ +unsigned char *scratch, *buffer; + + /* + * log2n is log2 of max(nqx,nqy) rounded up to next power of 2 + */ + nqmax = (nqx>nqy) ? nqx : nqy; + log2n = (int) (log((float) nqmax)/log(2.0)+0.5); + if (nqmax > (1<= 0; bit--) { + /* + * initial bit buffer + */ + b = 0; + bitbuffer = 0; + bits_to_go3 = 0; + /* + * on first pass copy A to scratch array + */ + qtree_onebit64(a,n,nqx,nqy,scratch,bit); + nx = (nqx+1)>>1; + ny = (nqy+1)>>1; + /* + * copy non-zero values to output buffer, which will be written + * in reverse order + */ + if (bufcopy(scratch,nx*ny,buffer,&b,bmax)) { + /* + * quadtree is expanding data, + * change warning code and just fill buffer with bit-map + */ + write_bdirect64(outfile,a,n,nqx,nqy,scratch,bit); + goto bitplane_done; + } + /* + * do log2n reductions + */ + for (k = 1; k>1; + ny = (ny+1)>>1; + if (bufcopy(scratch,nx*ny,buffer,&b,bmax)) { + write_bdirect64(outfile,a,n,nqx,nqy,scratch,bit); + goto bitplane_done; + } + } + /* + * OK, we've got the code in buffer + * Write quadtree warning code, then write buffer in reverse order + */ + output_nybble(outfile,0xF); + if (b==0) { + if (bits_to_go3>0) { + /* + * put out the last few bits + */ + output_nbits(outfile, bitbuffer & ((1<0) { + /* + * put out the last few bits + */ + output_nbits(outfile, bitbuffer & ((1<=0; i--) { + output_nbits(outfile,buffer[i],8); + } + } + bitplane_done: ; + } + free(buffer); + free(scratch); + return(0); +} + +/* ######################################################################### */ +/* + * copy non-zero codes from array to buffer + */ +static int +bufcopy(unsigned char a[], int n, unsigned char buffer[], int *b, int bmax) +{ +int i; + + for (i = 0; i < n; i++) { + if (a[i] != 0) { + /* + * add Huffman code for a[i] to buffer + */ + bitbuffer |= code[a[i]] << bits_to_go3; + bits_to_go3 += ncode[a[i]]; + if (bits_to_go3 >= 8) { + buffer[*b] = bitbuffer & 0xFF; + *b += 1; + /* + * return warning code if we fill buffer + */ + if (*b >= bmax) return(1); + bitbuffer >>= 8; + bits_to_go3 -= 8; + } + } + } + return(0); +} + +/* ######################################################################### */ +/* + * Do first quadtree reduction step on bit BIT of array A. + * Results put into B. + * + */ +static void +qtree_onebit(int a[], int n, int nx, int ny, unsigned char b[], int bit) +{ +int i, j, k; +int b0, b1, b2, b3; +int s10, s00; + + /* + * use selected bit to get amount to shift + */ + b0 = 1<> bit; + k += 1; + s00 += 2; + s10 += 2; + } + if (j < ny) { + /* + * row size is odd, do last element in row + * s00+1,s10+1 are off edge + */ + b[k] = ( ((a[s10 ]<<1) & b1) + | ((a[s00 ]<<3) & b3) ) >> bit; + k += 1; + } + } + if (i < nx) { + /* + * column size is odd, do last row + * s10,s10+1 are off edge + */ + s00 = n*i; + for (j = 0; j> bit; + k += 1; + s00 += 2; + } + if (j < ny) { + /* + * both row and column size are odd, do corner element + * s00+1, s10, s10+1 are off edge + */ + b[k] = ( ((a[s00 ]<<3) & b3) ) >> bit; + k += 1; + } + } +} +/* ######################################################################### */ +/* + * Do first quadtree reduction step on bit BIT of array A. + * Results put into B. + * + */ +static void +qtree_onebit64(LONGLONG a[], int n, int nx, int ny, unsigned char b[], int bit) +{ +int i, j, k; +LONGLONG b0, b1, b2, b3; +int s10, s00; + + /* + * use selected bit to get amount to shift + */ + b0 = ((LONGLONG) 1)<> bit); + k += 1; + s00 += 2; + s10 += 2; + } + if (j < ny) { + /* + * row size is odd, do last element in row + * s00+1,s10+1 are off edge + */ + b[k] = (unsigned char) (( ((a[s10 ]<<1) & b1) + | ((a[s00 ]<<3) & b3) ) >> bit); + k += 1; + } + } + if (i < nx) { + /* + * column size is odd, do last row + * s10,s10+1 are off edge + */ + s00 = n*i; + for (j = 0; j> bit); + k += 1; + s00 += 2; + } + if (j < ny) { + /* + * both row and column size are odd, do corner element + * s00+1, s10, s10+1 are off edge + */ + b[k] = (unsigned char) (( ((a[s00 ]<<3) & b3) ) >> bit); + k += 1; + } + } +} + +/* ######################################################################### */ +/* + * do one quadtree reduction step on array a + * results put into b (which may be the same as a) + */ +static void +qtree_reduce(unsigned char a[], int n, int nx, int ny, unsigned char b[]) +{ +int i, j, k; +int s10, s00; + + k = 0; /* k is index of b[i/2,j/2] */ + for (i = 0; i +#include +#include +#include + +// EAM : dropping fitsio prototypes +# define LONGLONG long long +# define DATA_COMPRESSION_ERR 413 /* error in imcompress routines */ +# define DATA_DECOMPRESSION_ERR 414 /* error in imcompress routines */ + +/* WDP added test to see if min and max are already defined */ +#ifndef min +#define min(a,b) (((a)<(b))?(a):(b)) +#endif +#ifndef max +#define max(a,b) (((a)>(b))?(a):(b)) +#endif + +#define input_nybble(infile) input_nbits(infile,4) + +static long nextchar; + +static int decode(unsigned char *infile, int *a, int *nx, int *ny, int *scale); +static int decode64(unsigned char *infile, LONGLONG *a, int *nx, int *ny, int *scale); +static int hinv(int a[], int nx, int ny, int smooth ,int scale); +static int hinv64(LONGLONG a[], int nx, int ny, int smooth ,int scale); +static void undigitize(int a[], int nx, int ny, int scale); +static void undigitize64(LONGLONG a[], int nx, int ny, int scale); +static void unshuffle(int a[], int n, int n2, int tmp[]); +static void unshuffle64(LONGLONG a[], int n, int n2, LONGLONG tmp[]); +static void hsmooth(int a[], int nxtop, int nytop, int ny, int scale); +static void hsmooth64(LONGLONG a[], int nxtop, int nytop, int ny, int scale); +static void qread(unsigned char *infile,char *a, int n); +static int readint(unsigned char *infile); +static LONGLONG readlonglong(unsigned char *infile); +static int dodecode(unsigned char *infile, int a[], int nx, int ny, unsigned char nbitplanes[3]); +static int dodecode64(unsigned char *infile, LONGLONG a[], int nx, int ny, unsigned char nbitplanes[3]); +static int qtree_decode(unsigned char *infile, int a[], int n, int nqx, int nqy, int nbitplanes); +static int qtree_decode64(unsigned char *infile, LONGLONG a[], int n, int nqx, int nqy, int nbitplanes); +static void start_inputing_bits(); +static int input_bit(unsigned char *infile); +static int input_nbits(unsigned char *infile, int n); + +static void qtree_expand(unsigned char *infile, unsigned char a[], int nx, int ny, unsigned char b[]); +static void qtree_bitins(unsigned char a[], int nx, int ny, int b[], int n, int bit); +static void qtree_bitins64(unsigned char a[], int nx, int ny, LONGLONG b[], int n, int bit); +static void qtree_copy(unsigned char a[], int nx, int ny, unsigned char b[], int n); +static void read_bdirect(unsigned char *infile, int a[], int n, int nqx, int nqy, unsigned char scratch[], int bit); +static void read_bdirect64(unsigned char *infile, LONGLONG a[], int n, int nqx, int nqy, unsigned char scratch[], int bit); +static int input_huffman(unsigned char *infile); +static int myread(unsigned char *file, char buffer[], int n); + +/* ---------------------------------------------------------------------- */ +int fits_hdecompress(unsigned char *input, int smooth, int *a, int *ny, int *nx, + int *scale, int *status) +{ + /* + decompress the input byte stream using the H-compress algorithm + + input - input array of compressed bytes + a - pre-allocated array to hold the output uncompressed image + nx - returned X axis size + ny - returned Y axis size + + NOTE: the nx and ny dimensions as defined within this code are reversed from + the usual FITS notation. ny is the fastest varying dimension, which is + usually considered the X axis in the FITS image display + + */ +int stat; + + if (*status > 0) return(*status); + + /* decode the input array */ + + stat = decode(input, a, nx, ny, scale); + *status = stat; + if (stat) return(*status); + + /* + * Un-Digitize + */ + undigitize(a, *nx, *ny, *scale); + + /* + * Inverse H-transform + */ + stat = hinv(a, *nx, *ny, smooth, *scale); + *status = stat; + + return(*status); +} +/* ---------------------------------------------------------------------- */ +int fits_hdecompress64(unsigned char *input, int smooth, LONGLONG *a, int *ny, int *nx, + int *scale, int *status) +{ + /* + decompress the input byte stream using the H-compress algorithm + + input - input array of compressed bytes + a - pre-allocated array to hold the output uncompressed image + nx - returned X axis size + ny - returned Y axis size + + NOTE: the nx and ny dimensions as defined within this code are reversed from + the usual FITS notation. ny is the fastest varying dimension, which is + usually considered the X axis in the FITS image display + + */ + int stat, *iarray, ii, nval; + + if (*status > 0) return(*status); + + /* decode the input array */ + + stat = decode64(input, a, nx, ny, scale); + *status = stat; + if (stat) return(*status); + + /* + * Un-Digitize + */ + undigitize64(a, *nx, *ny, *scale); + + /* + * Inverse H-transform + */ + stat = hinv64(a, *nx, *ny, smooth, *scale); + + *status = stat; + + /* pack the I*8 values back into an I*4 array */ + iarray = (int *) a; + nval = (*nx) * (*ny); + + for (ii = 0; ii < nval; ii++) + iarray[ii] = (int) a[ii]; + + return(*status); +} + +/* ############################################################################ */ +/* ############################################################################ */ + +/* Copyright (c) 1993 Association of Universities for Research + * in Astronomy. All rights reserved. Produced under National + * Aeronautics and Space Administration Contract No. NAS5-26555. + */ +/* hinv.c Inverse H-transform of NX x NY integer image + * + * Programmer: R. White Date: 23 July 1993 + */ + +/* ############################################################################ */ +static int +hinv(int a[], int nx, int ny, int smooth ,int scale) +/* +int smooth; 0 for no smoothing, else smooth during inversion +int scale; used if smoothing is specified +*/ +{ +int nmax, log2n, i, j, k; +int nxtop,nytop,nxf,nyf,c; +int oddx,oddy; +int shift, bit0, bit1, bit2, mask0, mask1, mask2, + prnd0, prnd1, prnd2, nrnd0, nrnd1, nrnd2, lowbit0, lowbit1; +int h0, hx, hy, hc; +int s10, s00; +int *tmp; + + /* + * log2n is log2 of max(nx,ny) rounded up to next power of 2 + */ + nmax = (nx>ny) ? nx : ny; + log2n = (int) (log((float) nmax)/log(2.0)+0.5); + if ( nmax > (1<> 1; + prnd1 = bit1 >> 1; + prnd2 = bit2 >> 1; + nrnd0 = prnd0 - 1; + nrnd1 = prnd1 - 1; + nrnd2 = prnd2 - 1; + /* + * round h0 to multiple of bit2 + */ + a[0] = (a[0] + ((a[0] >= 0) ? prnd2 : nrnd2)) & mask2; + /* + * do log2n expansions + * + * We're indexing a as a 2-D array with dimensions (nx,ny). + */ + nxtop = 1; + nytop = 1; + nxf = nx; + nyf = ny; + c = 1<=0; k--) { + /* + * this somewhat cryptic code generates the sequence + * ntop[k-1] = (ntop[k]+1)/2, where ntop[log2n] = n + */ + c = c>>1; + nxtop = nxtop<<1; + nytop = nytop<<1; + if (nxf <= c) { nxtop -= 1; } else { nxf -= c; } + if (nyf <= c) { nytop -= 1; } else { nyf -= c; } + /* + * double shift and fix nrnd0 (because prnd0=0) on last pass + */ + if (k == 0) { + nrnd0 = 0; + shift = 2; + } + /* + * unshuffle in each dimension to interleave coefficients + */ + for (i = 0; i= 0) ? prnd1 : nrnd1)) & mask1; + hy = (hy + ((hy >= 0) ? prnd1 : nrnd1)) & mask1; + hc = (hc + ((hc >= 0) ? prnd0 : nrnd0)) & mask0; + /* + * propagate bit0 of hc to hx,hy + */ + lowbit0 = hc & bit0; + hx = (hx >= 0) ? (hx - lowbit0) : (hx + lowbit0); + hy = (hy >= 0) ? (hy - lowbit0) : (hy + lowbit0); + /* + * Propagate bits 0 and 1 of hc,hx,hy to h0. + * This could be simplified if we assume h0>0, but then + * the inversion would not be lossless for images with + * negative pixels. + */ + lowbit1 = (hc ^ hx ^ hy) & bit1; + h0 = (h0 >= 0) + ? (h0 + lowbit0 - lowbit1) + : (h0 + ((lowbit0 == 0) ? lowbit1 : (lowbit0-lowbit1))); + /* + * Divide sums by 2 (4 last time) + */ + a[s10+1] = (h0 + hx + hy + hc) >> shift; + a[s10 ] = (h0 + hx - hy - hc) >> shift; + a[s00+1] = (h0 - hx + hy - hc) >> shift; + a[s00 ] = (h0 - hx - hy + hc) >> shift; + s00 += 2; + s10 += 2; + } + if (oddy) { + /* + * do last element in row if row length is odd + * s00+1, s10+1 are off edge + */ + h0 = a[s00 ]; + hx = a[s10 ]; + hx = ((hx >= 0) ? (hx+prnd1) : (hx+nrnd1)) & mask1; + lowbit1 = hx & bit1; + h0 = (h0 >= 0) ? (h0 - lowbit1) : (h0 + lowbit1); + a[s10 ] = (h0 + hx) >> shift; + a[s00 ] = (h0 - hx) >> shift; + } + } + if (oddx) { + /* + * do last row if column length is odd + * s10, s10+1 are off edge + */ + s00 = ny*i; + for (j = 0; j= 0) ? (hy+prnd1) : (hy+nrnd1)) & mask1; + lowbit1 = hy & bit1; + h0 = (h0 >= 0) ? (h0 - lowbit1) : (h0 + lowbit1); + a[s00+1] = (h0 + hy) >> shift; + a[s00 ] = (h0 - hy) >> shift; + s00 += 2; + } + if (oddy) { + /* + * do corner element if both row and column lengths are odd + * s00+1, s10, s10+1 are off edge + */ + h0 = a[s00 ]; + a[s00 ] = h0 >> shift; + } + } + /* + * divide all the masks and rounding values by 2 + */ + bit2 = bit1; + bit1 = bit0; + bit0 = bit0 >> 1; + mask1 = mask0; + mask0 = mask0 >> 1; + prnd1 = prnd0; + prnd0 = prnd0 >> 1; + nrnd1 = nrnd0; + nrnd0 = prnd0 - 1; + } + free(tmp); + return(0); +} +/* ############################################################################ */ +static int +hinv64(LONGLONG a[], int nx, int ny, int smooth ,int scale) +/* +int smooth; 0 for no smoothing, else smooth during inversion +int scale; used if smoothing is specified +*/ +{ +int nmax, log2n, i, j, k; +int nxtop,nytop,nxf,nyf,c; +int oddx,oddy; +int shift; +LONGLONG mask0, mask1, mask2, prnd0, prnd1, prnd2, bit0, bit1, bit2; +LONGLONG nrnd0, nrnd1, nrnd2, lowbit0, lowbit1; +LONGLONG h0, hx, hy, hc; +int s10, s00; +LONGLONG *tmp; + + /* + * log2n is log2 of max(nx,ny) rounded up to next power of 2 + */ + nmax = (nx>ny) ? nx : ny; + log2n = (int) (log((float) nmax)/log(2.0)+0.5); + if ( nmax > (1<> 1; + prnd1 = bit1 >> 1; + prnd2 = bit2 >> 1; + nrnd0 = prnd0 - 1; + nrnd1 = prnd1 - 1; + nrnd2 = prnd2 - 1; + /* + * round h0 to multiple of bit2 + */ + a[0] = (a[0] + ((a[0] >= 0) ? prnd2 : nrnd2)) & mask2; + /* + * do log2n expansions + * + * We're indexing a as a 2-D array with dimensions (nx,ny). + */ + nxtop = 1; + nytop = 1; + nxf = nx; + nyf = ny; + c = 1<=0; k--) { + /* + * this somewhat cryptic code generates the sequence + * ntop[k-1] = (ntop[k]+1)/2, where ntop[log2n] = n + */ + c = c>>1; + nxtop = nxtop<<1; + nytop = nytop<<1; + if (nxf <= c) { nxtop -= 1; } else { nxf -= c; } + if (nyf <= c) { nytop -= 1; } else { nyf -= c; } + /* + * double shift and fix nrnd0 (because prnd0=0) on last pass + */ + if (k == 0) { + nrnd0 = 0; + shift = 2; + } + /* + * unshuffle in each dimension to interleave coefficients + */ + for (i = 0; i= 0) ? prnd1 : nrnd1)) & mask1; + hy = (hy + ((hy >= 0) ? prnd1 : nrnd1)) & mask1; + hc = (hc + ((hc >= 0) ? prnd0 : nrnd0)) & mask0; + /* + * propagate bit0 of hc to hx,hy + */ + lowbit0 = hc & bit0; + hx = (hx >= 0) ? (hx - lowbit0) : (hx + lowbit0); + hy = (hy >= 0) ? (hy - lowbit0) : (hy + lowbit0); + /* + * Propagate bits 0 and 1 of hc,hx,hy to h0. + * This could be simplified if we assume h0>0, but then + * the inversion would not be lossless for images with + * negative pixels. + */ + lowbit1 = (hc ^ hx ^ hy) & bit1; + h0 = (h0 >= 0) + ? (h0 + lowbit0 - lowbit1) + : (h0 + ((lowbit0 == 0) ? lowbit1 : (lowbit0-lowbit1))); + /* + * Divide sums by 2 (4 last time) + */ + a[s10+1] = (h0 + hx + hy + hc) >> shift; + a[s10 ] = (h0 + hx - hy - hc) >> shift; + a[s00+1] = (h0 - hx + hy - hc) >> shift; + a[s00 ] = (h0 - hx - hy + hc) >> shift; + s00 += 2; + s10 += 2; + } + if (oddy) { + /* + * do last element in row if row length is odd + * s00+1, s10+1 are off edge + */ + h0 = a[s00 ]; + hx = a[s10 ]; + hx = ((hx >= 0) ? (hx+prnd1) : (hx+nrnd1)) & mask1; + lowbit1 = hx & bit1; + h0 = (h0 >= 0) ? (h0 - lowbit1) : (h0 + lowbit1); + a[s10 ] = (h0 + hx) >> shift; + a[s00 ] = (h0 - hx) >> shift; + } + } + if (oddx) { + /* + * do last row if column length is odd + * s10, s10+1 are off edge + */ + s00 = ny*i; + for (j = 0; j= 0) ? (hy+prnd1) : (hy+nrnd1)) & mask1; + lowbit1 = hy & bit1; + h0 = (h0 >= 0) ? (h0 - lowbit1) : (h0 + lowbit1); + a[s00+1] = (h0 + hy) >> shift; + a[s00 ] = (h0 - hy) >> shift; + s00 += 2; + } + if (oddy) { + /* + * do corner element if both row and column lengths are odd + * s00+1, s10, s10+1 are off edge + */ + h0 = a[s00 ]; + a[s00 ] = h0 >> shift; + } + } + /* + * divide all the masks and rounding values by 2 + */ + bit2 = bit1; + bit1 = bit0; + bit0 = bit0 >> 1; + mask1 = mask0; + mask0 = mask0 >> 1; + prnd1 = prnd0; + prnd0 = prnd0 >> 1; + nrnd1 = nrnd0; + nrnd0 = prnd0 - 1; + } + free(tmp); + return(0); +} + +/* ############################################################################ */ +static void +unshuffle(int a[], int n, int n2, int tmp[]) +/* +int a[]; array to shuffle +int n; number of elements to shuffle +int n2; second dimension +int tmp[]; scratch storage +*/ +{ +int i; +int nhalf; +int *p1, *p2, *pt; + + /* + * copy 2nd half of array to tmp + */ + nhalf = (n+1)>>1; + pt = tmp; + p1 = &a[n2*nhalf]; /* pointer to a[i] */ + for (i=nhalf; i= 0; i--) { + *p1 = *p2; + p2 -= n2; + p1 -= (n2+n2); + } + /* + * now distribute 2nd half of array (in tmp) to odd elements + */ + pt = tmp; + p1 = &a[n2]; /* pointer to a[i] */ + for (i=1; i>1; + pt = tmp; + p1 = &a[n2*nhalf]; /* pointer to a[i] */ + for (i=nhalf; i= 0; i--) { + *p1 = *p2; + p2 -= n2; + p1 -= (n2+n2); + } + /* + * now distribute 2nd half of array (in tmp) to odd elements + */ + pt = tmp; + p1 = &a[n2]; /* pointer to a[i] */ + for (i=1; i> 1); + if (smax <= 0) return; + ny2 = ny << 1; + /* + * We're indexing a as a 2-D array with dimensions (nxtop,ny) of which + * only (nxtop,nytop) are used. The coefficients on the edge of the + * array are not adjusted (which is why the loops below start at 2 + * instead of 0 and end at nxtop-2 instead of nxtop.) + */ + /* + * Adjust x difference hx + */ + for (i = 2; i=0, dmin<=0. + */ + if (dmin < dmax) { + diff = max( min(diff, dmax), dmin); + /* + * Compute change in slope limited to range +/- smax. + * Careful with rounding negative numbers when using + * shift for divide by 8. + */ + s = diff-(a[s10]<<3); + s = (s>=0) ? (s>>3) : ((s+7)>>3) ; + s = max( min(s, smax), -smax); + a[s10] = a[s10]+s; + } + s00 += 2; + s10 += 2; + } + } + /* + * Adjust y difference hy + */ + for (i = 0; i=0) ? (s>>3) : ((s+7)>>3) ; + s = max( min(s, smax), -smax); + a[s00+1] = a[s00+1]+s; + } + s00 += 2; + s10 += 2; + } + } + /* + * Adjust curvature difference hc + */ + for (i = 2; i=0, dmin<=0. + */ + if (dmin < dmax) { + diff = max( min(diff, dmax), dmin); + /* + * Compute change in slope limited to range +/- smax. + * Careful with rounding negative numbers when using + * shift for divide by 64. + */ + s = diff-(a[s10+1]<<6); + s = (s>=0) ? (s>>6) : ((s+63)>>6) ; + s = max( min(s, smax), -smax); + a[s10+1] = a[s10+1]+s; + } + s00 += 2; + s10 += 2; + } + } +} +/* ############################################################################ */ +static void +hsmooth64(LONGLONG a[], int nxtop, int nytop, int ny, int scale) +/* +LONGLONG a[]; array of H-transform coefficients +int nxtop,nytop; size of coefficient block to use +int ny; actual 1st dimension of array +int scale; truncation scale factor that was used +*/ +{ +int i, j; +int ny2, s10, s00; +LONGLONG hm, h0, hp, hmm, hpm, hmp, hpp, hx2, hy2, diff, dmax, dmin, s, smax, m1, m2; + + /* + * Maximum change in coefficients is determined by scale factor. + * Since we rounded during division (see digitize.c), the biggest + * permitted change is scale/2. + */ + smax = (scale >> 1); + if (smax <= 0) return; + ny2 = ny << 1; + /* + * We're indexing a as a 2-D array with dimensions (nxtop,ny) of which + * only (nxtop,nytop) are used. The coefficients on the edge of the + * array are not adjusted (which is why the loops below start at 2 + * instead of 0 and end at nxtop-2 instead of nxtop.) + */ + /* + * Adjust x difference hx + */ + for (i = 2; i=0, dmin<=0. + */ + if (dmin < dmax) { + diff = max( min(diff, dmax), dmin); + /* + * Compute change in slope limited to range +/- smax. + * Careful with rounding negative numbers when using + * shift for divide by 8. + */ + s = diff-(a[s10]<<3); + s = (s>=0) ? (s>>3) : ((s+7)>>3) ; + s = max( min(s, smax), -smax); + a[s10] = a[s10]+s; + } + s00 += 2; + s10 += 2; + } + } + /* + * Adjust y difference hy + */ + for (i = 0; i=0) ? (s>>3) : ((s+7)>>3) ; + s = max( min(s, smax), -smax); + a[s00+1] = a[s00+1]+s; + } + s00 += 2; + s10 += 2; + } + } + /* + * Adjust curvature difference hc + */ + for (i = 2; i=0, dmin<=0. + */ + if (dmin < dmax) { + diff = max( min(diff, dmax), dmin); + /* + * Compute change in slope limited to range +/- smax. + * Careful with rounding negative numbers when using + * shift for divide by 64. + */ + s = diff-(a[s10+1]<<6); + s = (s>=0) ? (s>>6) : ((s+63)>>6) ; + s = max( min(s, smax), -smax); + a[s10+1] = a[s10+1]+s; + } + s00 += 2; + s10 += 2; + } + } +} + + +/* ############################################################################ */ +/* ############################################################################ */ +/* Copyright (c) 1993 Association of Universities for Research + * in Astronomy. All rights reserved. Produced under National + * Aeronautics and Space Administration Contract No. NAS5-26555. + */ +/* undigitize.c undigitize H-transform + * + * Programmer: R. White Date: 9 May 1991 + */ + +/* ############################################################################ */ +static void +undigitize(int a[], int nx, int ny, int scale) +{ +int *p; + + /* + * multiply by scale + */ + if (scale <= 1) return; + for (p=a; p <= &a[nx*ny-1]; p++) *p = (*p)*scale; +} +/* ############################################################################ */ +static void +undigitize64(LONGLONG a[], int nx, int ny, int scale) +{ +LONGLONG *p, scale64; + + /* + * multiply by scale + */ + if (scale <= 1) return; + scale64 = (LONGLONG) scale; /* use a 64-bit int for efficiency in the big loop */ + + for (p=a; p <= &a[nx*ny-1]; p++) *p = (*p)*scale64; +} + +/* ############################################################################ */ +/* ############################################################################ */ +/* Copyright (c) 1993 Association of Universities for Research + * in Astronomy. All rights reserved. Produced under National + * Aeronautics and Space Administration Contract No. NAS5-26555. + */ +/* decode.c read codes from infile and construct array + * + * Programmer: R. White Date: 2 February 1994 + */ + + +static char code_magic[2] = { (char)0xDD, (char)0x99 }; + +/* ############################################################################ */ +static int decode(unsigned char *infile, int *a, int *nx, int *ny, int *scale) +/* +char *infile; input file +int *a; address of output array [nx][ny] +int *nx,*ny; size of output array +int *scale; scale factor for digitization +*/ +{ +LONGLONG sumall; +int nel, stat; +unsigned char nbitplanes[3]; +char tmagic[2]; + + /* initialize the byte read position to the beginning of the array */; + nextchar = 0; + + /* + * File starts either with special 2-byte magic code or with + * FITS keyword "SIMPLE =" + */ + qread(infile, tmagic, sizeof(tmagic)); + /* + * check for correct magic code value + */ + if (memcmp(tmagic,code_magic,sizeof(code_magic)) != 0) { + fprintf (stderr, "bad file format"); + return(DATA_DECOMPRESSION_ERR); + } + *nx =readint(infile); /* x size of image */ + *ny =readint(infile); /* y size of image */ + *scale=readint(infile); /* scale factor for digitization */ + + nel = (*nx) * (*ny); + + /* sum of all pixels */ + sumall=readlonglong(infile); + /* # bits in quadrants */ + + qread(infile, (char *) nbitplanes, sizeof(nbitplanes)); + + stat = dodecode(infile, a, *nx, *ny, nbitplanes); + /* + * put sum of all pixels back into pixel 0 + */ + a[0] = (int) sumall; + return(stat); +} +/* ############################################################################ */ +static int decode64(unsigned char *infile, LONGLONG *a, int *nx, int *ny, int *scale) +/* +char *infile; input file +LONGLONG *a; address of output array [nx][ny] +int *nx,*ny; size of output array +int *scale; scale factor for digitization +*/ +{ +int nel, stat; +LONGLONG sumall; +unsigned char nbitplanes[3]; +char tmagic[2]; + + /* initialize the byte read position to the beginning of the array */; + nextchar = 0; + + /* + * File starts either with special 2-byte magic code or with + * FITS keyword "SIMPLE =" + */ + qread(infile, tmagic, sizeof(tmagic)); + /* + * check for correct magic code value + */ + if (memcmp(tmagic,code_magic,sizeof(code_magic)) != 0) { + fprintf (stderr, "bad file format"); + return(DATA_DECOMPRESSION_ERR); + } + *nx =readint(infile); /* x size of image */ + *ny =readint(infile); /* y size of image */ + *scale=readint(infile); /* scale factor for digitization */ + + nel = (*nx) * (*ny); + + /* sum of all pixels */ + sumall=readlonglong(infile); + /* # bits in quadrants */ + + qread(infile, (char *) nbitplanes, sizeof(nbitplanes)); + + stat = dodecode64(infile, a, *nx, *ny, nbitplanes); + /* + * put sum of all pixels back into pixel 0 + */ + a[0] = sumall; + + return(stat); +} + + +/* ############################################################################ */ +/* ############################################################################ */ +/* Copyright (c) 1993 Association of Universities for Research + * in Astronomy. All rights reserved. Produced under National + * Aeronautics and Space Administration Contract No. NAS5-26555. + */ +/* dodecode.c Decode stream of characters on infile and return array + * + * This version encodes the different quadrants separately + * + * Programmer: R. White Date: 9 May 1991 + */ + +/* ############################################################################ */ +static int +dodecode(unsigned char *infile, int a[], int nx, int ny, unsigned char nbitplanes[3]) + +/* int a[]; + int nx,ny; Array dimensions are [nx][ny] + unsigned char nbitplanes[3]; Number of bit planes in quadrants +*/ +{ +int i, nel, nx2, ny2, stat; + + nel = nx*ny; + nx2 = (nx+1)/2; + ny2 = (ny+1)/2; + + /* + * initialize a to zero + */ + for (i=0; inqy) ? nqx : nqy; + log2n = (int) (log((float) nqmax)/log(2.0)+0.5); + if (nqmax > (1<= 0; bit--) { + /* + * Was bitplane was quadtree-coded or written directly? + */ + b = input_nybble(infile); + + if(b == 0) { + /* + * bit map was written directly + */ + read_bdirect(infile,a,n,nqx,nqy,scratch,bit); + } else if (b != 0xf) { + fprintf (stderr, "qtree_decode: bad format code"); + return(DATA_DECOMPRESSION_ERR); + } else { + /* + * bitmap was quadtree-coded, do log2n expansions + * + * read first code + */ + scratch[0] = input_huffman(infile); + /* + * now do log2n expansions, reading codes from file as necessary + */ + nx = 1; + ny = 1; + nfx = nqx; + nfy = nqy; + c = 1<>1; + nx = nx<<1; + ny = ny<<1; + if (nfx <= c) { nx -= 1; } else { nfx -= c; } + if (nfy <= c) { ny -= 1; } else { nfy -= c; } + qtree_expand(infile,scratch,nx,ny,scratch); + } + /* + * now copy last set of 4-bit codes to bitplane bit of array a + */ + qtree_bitins(scratch,nqx,nqy,a,n,bit); + } + } + free(scratch); + return(0); +} +/* ############################################################################ */ +static int +qtree_decode64(unsigned char *infile, LONGLONG a[], int n, int nqx, int nqy, int nbitplanes) + +/* +char *infile; +LONGLONG a[]; a is 2-D array with dimensions (n,n) +int n; length of full row in a +int nqx; partial length of row to decode +int nqy; partial length of column (<=n) +int nbitplanes; number of bitplanes to decode +*/ +{ +int log2n, k, bit, b, nqmax; +int nx,ny,nfx,nfy,c; +int nqx2, nqy2; +unsigned char *scratch; + + /* + * log2n is log2 of max(nqx,nqy) rounded up to next power of 2 + */ + nqmax = (nqx>nqy) ? nqx : nqy; + log2n = (int) (log((float) nqmax)/log(2.0)+0.5); + if (nqmax > (1<= 0; bit--) { + /* + * Was bitplane was quadtree-coded or written directly? + */ + b = input_nybble(infile); + + if(b == 0) { + /* + * bit map was written directly + */ + read_bdirect64(infile,a,n,nqx,nqy,scratch,bit); + } else if (b != 0xf) { + fprintf (stderr, "qtree_decode64: bad format code"); + return(DATA_DECOMPRESSION_ERR); + } else { + /* + * bitmap was quadtree-coded, do log2n expansions + * + * read first code + */ + scratch[0] = input_huffman(infile); + /* + * now do log2n expansions, reading codes from file as necessary + */ + nx = 1; + ny = 1; + nfx = nqx; + nfy = nqy; + c = 1<>1; + nx = nx<<1; + ny = ny<<1; + if (nfx <= c) { nx -= 1; } else { nfx -= c; } + if (nfy <= c) { ny -= 1; } else { nfy -= c; } + qtree_expand(infile,scratch,nx,ny,scratch); + } + /* + * now copy last set of 4-bit codes to bitplane bit of array a + */ + qtree_bitins64(scratch,nqx,nqy,a,n,bit); + } + } + free(scratch); + return(0); +} + + +/* ############################################################################ */ +/* + * do one quadtree expansion step on array a[(nqx+1)/2,(nqy+1)/2] + * results put into b[nqx,nqy] (which may be the same as a) + */ +static void +qtree_expand(unsigned char *infile, unsigned char a[], int nx, int ny, unsigned char b[]) +{ +int i; + + /* + * first copy a to b, expanding each 4-bit value + */ + qtree_copy(a,nx,ny,b,ny); + /* + * now read new 4-bit values into b for each non-zero element + */ + for (i = nx*ny-1; i >= 0; i--) { + if (b[i] != 0) b[i] = input_huffman(infile); + } +} + +/* ############################################################################ */ +/* + * copy 4-bit values from a[(nx+1)/2,(ny+1)/2] to b[nx,ny], expanding + * each value to 2x2 pixels + * a,b may be same array + */ +static void +qtree_copy(unsigned char a[], int nx, int ny, unsigned char b[], int n) +/* int n; declared y dimension of b */ +{ +int i, j, k, nx2, ny2; +int s00, s10; + + /* + * first copy 4-bit values to b + * start at end in case a,b are same array + */ + nx2 = (nx+1)/2; + ny2 = (ny+1)/2; + k = ny2*(nx2-1)+ny2-1; /* k is index of a[i,j] */ + for (i = nx2-1; i >= 0; i--) { + s00 = 2*(n*i+ny2-1); /* s00 is index of b[2*i,2*j] */ + for (j = ny2-1; j >= 0; j--) { + b[s00] = a[k]; + k -= 1; + s00 -= 2; + } + } + /* + * now expand each 2x2 block + */ + for (i = 0; i>1) & 1; + b[s00+1] = (b[s00]>>2) & 1; + b[s00 ] = (b[s00]>>3) & 1; + s00 += 2; + s10 += 2; + } + if (j < ny) { + /* + * row size is odd, do last element in row + * s00+1, s10+1 are off edge + */ + b[s10 ] = (b[s00]>>1) & 1; + b[s00 ] = (b[s00]>>3) & 1; + } + } + if (i < nx) { + /* + * column size is odd, do last row + * s10, s10+1 are off edge + */ + s00 = n*i; + for (j = 0; j>2) & 1; + b[s00 ] = (b[s00]>>3) & 1; + s00 += 2; + } + if (j < ny) { + /* + * both row and column size are odd, do corner element + * s00+1, s10, s10+1 are off edge + */ + b[s00 ] = (b[s00]>>3) & 1; + } + } +} + +/* ############################################################################ */ +/* + * Copy 4-bit values from a[(nx+1)/2,(ny+1)/2] to b[nx,ny], expanding + * each value to 2x2 pixels and inserting into bitplane BIT of B. + * A,B may NOT be same array (it wouldn't make sense to be inserting + * bits into the same array anyway.) + */ +static void +qtree_bitins(unsigned char a[], int nx, int ny, int b[], int n, int bit) +/* + int n; declared y dimension of b +*/ +{ +int i, j, k; +int s00, s10; + + /* + * expand each 2x2 block + */ + k = 0; /* k is index of a[i/2,j/2] */ + for (i = 0; i>1) & 1) << bit; + b[s00+1] |= ((a[k]>>2) & 1) << bit; + b[s00 ] |= ((a[k]>>3) & 1) << bit; + s00 += 2; + s10 += 2; + k += 1; + } + if (j < ny) { + /* + * row size is odd, do last element in row + * s00+1, s10+1 are off edge + */ + b[s10 ] |= ((a[k]>>1) & 1) << bit; + b[s00 ] |= ((a[k]>>3) & 1) << bit; + k += 1; + } + } + if (i < nx) { + /* + * column size is odd, do last row + * s10, s10+1 are off edge + */ + s00 = n*i; + for (j = 0; j>2) & 1) << bit; + b[s00 ] |= ((a[k]>>3) & 1) << bit; + s00 += 2; + k += 1; + } + if (j < ny) { + /* + * both row and column size are odd, do corner element + * s00+1, s10, s10+1 are off edge + */ + b[s00 ] |= ((a[k]>>3) & 1) << bit; + k += 1; + } + } +} +/* ############################################################################ */ +/* + * Copy 4-bit values from a[(nx+1)/2,(ny+1)/2] to b[nx,ny], expanding + * each value to 2x2 pixels and inserting into bitplane BIT of B. + * A,B may NOT be same array (it wouldn't make sense to be inserting + * bits into the same array anyway.) + */ +static void +qtree_bitins64(unsigned char a[], int nx, int ny, LONGLONG b[], int n, int bit) +/* + int n; declared y dimension of b +*/ +{ +int i, j, k; +int s00, s10; + + /* + * expand each 2x2 block + */ + k = 0; /* k is index of a[i/2,j/2] */ + for (i = 0; i>1) & 1) << bit; + b[s00+1] |= ((((LONGLONG)a[k])>>2) & 1) << bit; + b[s00 ] |= ((((LONGLONG)a[k])>>3) & 1) << bit; + s00 += 2; + s10 += 2; + k += 1; + } + if (j < ny) { + /* + * row size is odd, do last element in row + * s00+1, s10+1 are off edge + */ + b[s10 ] |= ((((LONGLONG)a[k])>>1) & 1) << bit; + b[s00 ] |= ((((LONGLONG)a[k])>>3) & 1) << bit; + k += 1; + } + } + if (i < nx) { + /* + * column size is odd, do last row + * s10, s10+1 are off edge + */ + s00 = n*i; + for (j = 0; j>2) & 1) << bit; + b[s00 ] |= ((((LONGLONG)a[k])>>3) & 1) << bit; + s00 += 2; + k += 1; + } + if (j < ny) { + /* + * both row and column size are odd, do corner element + * s00+1, s10, s10+1 are off edge + */ + b[s00 ] |= ((((LONGLONG)a[k])>>3) & 1) << bit; + k += 1; + } + } +} + +/* ############################################################################ */ +static void +read_bdirect(unsigned char *infile, int a[], int n, int nqx, int nqy, unsigned char scratch[], int bit) +{ +int i; + + /* + * read bit image packed 4 pixels/nybble + */ + for (i = 0; i < ((nqx+1)/2) * ((nqy+1)/2); i++) { + scratch[i] = input_nybble(infile); + } + /* + * insert in bitplane BIT of image A + */ + qtree_bitins(scratch,nqx,nqy,a,n,bit); +} +/* ############################################################################ */ +static void +read_bdirect64(unsigned char *infile, LONGLONG a[], int n, int nqx, int nqy, unsigned char scratch[], int bit) +{ +int i; + + /* + * read bit image packed 4 pixels/nybble + */ + for (i = 0; i < ((nqx+1)/2) * ((nqy+1)/2); i++) { + scratch[i] = input_nybble(infile); + } + /* + * insert in bitplane BIT of image A + */ + qtree_bitins64(scratch,nqx,nqy,a,n,bit); +} + +/* ############################################################################ */ +/* + * Huffman decoding for fixed codes + * + * Coded values range from 0-15 + * + * Huffman code values (hex): + * + * 3e, 00, 01, 08, 02, 09, 1a, 1b, + * 03, 1c, 0a, 1d, 0b, 1e, 3f, 0c + * + * and number of bits in each code: + * + * 6, 3, 3, 4, 3, 4, 5, 5, + * 3, 5, 4, 5, 4, 5, 6, 4 + */ +static int input_huffman(unsigned char *infile) +{ +int c; + + /* + * get first 3 bits to start + */ + c = input_nbits(infile,3); + if (c < 4) { + /* + * this is all we need + * return 1,2,4,8 for c=0,1,2,3 + */ + return(1< 0) { + total += nread; + if (total==n) return(total); + } +*/ + + memcpy(buffer, &file[nextchar], n); + nextchar += n; + + return(n); +} + +/* ############################################################################ */ +/* ############################################################################ */ +/* Copyright (c) 1993 Association of Universities for Research + * in Astronomy. All rights reserved. Produced under National + * Aeronautics and Space Administration Contract No. NAS5-26555. + */ + +/* BIT INPUT ROUTINES */ + +/* THE BIT BUFFER */ + +static int buffer2; /* Bits waiting to be input */ +static int bits_to_go; /* Number of bits still in buffer */ + + +/* INITIALIZE BIT INPUT */ + +/* ############################################################################ */ +static void start_inputing_bits() +{ + /* + * Buffer starts out with no bits in it + */ + bits_to_go = 0; +} + + +/* ############################################################################ */ +/* INPUT A BIT */ + +static int input_bit(unsigned char *infile) +{ + if (bits_to_go == 0) { /* Read the next byte if no */ + /* buffer2 = getc(infile); */ /* bits are left in buffer */ + + buffer2 = infile[nextchar]; + nextchar++; + + /* if (buffer2 == EOF) { */ + /* + * end of file is an error for this application + */ +/* + fprintf(stderr, "input_bit: unexpected end-of-file\n"); + exit(-1); + } +*/ + bits_to_go = 8; + } + /* + * Return the next bit + */ + bits_to_go -= 1; + return((buffer2>>bits_to_go) & 1); +} + + +/* ############################################################################ */ +/* INPUT N BITS (N must be <= 8) */ + +static int input_nbits(unsigned char *infile, int n) +{ +int c; + + if (bits_to_go < n) { + /* + * need another byte's worth of bits + */ + buffer2 <<= 8; + +/* c = getc(infile); */ + + c = infile[nextchar]; + nextchar++; + +/* if (c == EOF) { */ + /* + * end of file is an error for this application + */ +/* + fprintf(stderr, "input_nbits: unexpected end-of-file\n"); + exit(-1); + } +*/ + buffer2 |= c; + bits_to_go += 8; + } + /* + * now pick off the first n bits + */ + bits_to_go -= n; + return( (buffer2>>bits_to_go) & ((1< +/* + The following prototype code was provided by Doug Tody, NRAO, for + performing conversion between pixel arrays and line lists. The + compression technique is used in IRAF. +*/ +int pl_p2li (int *pxsrc, int xs, short *lldst, int npix); +int pl_l2pi (short *ll_src, int xs, int *px_dst, int npix); + + +/* + * PL_P2L -- Convert a pixel array to a line list. The length of the list is + * returned as the function value. + * + * Translated from the SPP version using xc -f, f2c. 8Sep99 DCT. + */ + +#ifndef min +#define min(a,b) (((a)<(b))?(a):(b)) +#endif +#ifndef max +#define max(a,b) (((a)>(b))?(a):(b)) +#endif + +int pl_p2li (int *pxsrc, int xs, short *lldst, int npix) +/* int *pxsrc; input pixel array */ +/* int xs; starting index in pxsrc (?) */ +/* short *lldst; encoded line list */ +/* int npix; number of pixels to convert */ +{ + /* System generated locals */ + int ret_val, i__1, i__2, i__3; + + /* Local variables */ + static int zero, v, x1, hi, ip, dv, xe, np, op, iz, nv, pv, nz; + + /* Parameter adjustments */ + --lldst; + --pxsrc; + + /* Function Body */ + if (! (npix <= 0)) { + goto L110; + } + ret_val = 0; + goto L100; +L110: + lldst[3] = -100; + lldst[2] = 7; + lldst[1] = 0; + lldst[6] = 0; + lldst[7] = 0; + xe = xs + npix - 1; + op = 8; + zero = 0; +/* Computing MAX */ + i__1 = zero, i__2 = pxsrc[xs]; + pv = max(i__1,i__2); + x1 = xs; + iz = xs; + hi = 1; + i__1 = xe; + for (ip = xs; ip <= i__1; ++ip) { + if (! (ip < xe)) { + goto L130; + } +/* Computing MAX */ + i__2 = zero, i__3 = pxsrc[ip + 1]; + nv = max(i__2,i__3); + if (! (nv == pv)) { + goto L140; + } + goto L120; +L140: + if (! (pv == 0)) { + goto L150; + } + pv = nv; + x1 = ip + 1; + goto L120; +L150: + goto L131; +L130: + if (! (pv == 0)) { + goto L160; + } + x1 = xe + 1; +L160: +L131: + np = ip - x1 + 1; + nz = x1 - iz; + if (! (pv > 0)) { + goto L170; + } + dv = pv - hi; + if (! (dv != 0)) { + goto L180; + } + hi = pv; + if (! (abs(dv) > 4095)) { + goto L190; + } + lldst[op] = (short) ((pv & 4095) + 4096); + ++op; + lldst[op] = (short) (pv / 4096); + ++op; + goto L191; +L190: + if (! (dv < 0)) { + goto L200; + } + lldst[op] = (short) (-dv + 12288); + goto L201; +L200: + lldst[op] = (short) (dv + 8192); +L201: + ++op; + if (! (np == 1 && nz == 0)) { + goto L210; + } + v = lldst[op - 1]; + lldst[op - 1] = (short) (v | 16384); + goto L91; +L210: +L191: +L180: +L170: + if (! (nz > 0)) { + goto L220; + } +L230: + if (! (nz > 0)) { + goto L232; + } + lldst[op] = (short) min(4095,nz); + ++op; +/* L231: */ + nz += -4095; + goto L230; +L232: + if (! (np == 1 && pv > 0)) { + goto L240; + } + lldst[op - 1] = (short) (lldst[op - 1] + 20481); + goto L91; +L240: +L220: +L250: + if (! (np > 0)) { + goto L252; + } + lldst[op] = (short) (min(4095,np) + 16384); + ++op; +/* L251: */ + np += -4095; + goto L250; +L252: +L91: + x1 = ip + 1; + iz = x1; + pv = nv; +L120: + ; + } +/* L121: */ + lldst[4] = (short) ((op - 1) % 32768); + lldst[5] = (short) ((op - 1) / 32768); + ret_val = op - 1; + goto L100; +L100: + return ret_val; +} /* plp2li_ */ + +/* + * PL_L2PI -- Translate a PLIO line list into an integer pixel array. + * The number of pixels output (always npix) is returned as the function + * value. + * + * Translated from the SPP version using xc -f, f2c. 8Sep99 DCT. + */ + +int pl_l2pi (short *ll_src, int xs, int *px_dst, int npix) +/* short *ll_src; encoded line list */ +/* int xs; starting index in ll_src */ +/* int *px_dst; output pixel array */ +/* int npix; number of pixels to convert */ +{ + /* System generated locals */ + int ret_val, i__1, i__2; + + /* Local variables */ + static int data, sw0001, otop, i__, lllen, i1, i2, x1, x2, ip, xe, np, + op, pv, opcode, llfirt; + static int skipwd; + + /* Parameter adjustments */ + --px_dst; + --ll_src; + + /* Function Body */ + if (! (ll_src[3] > 0)) { + goto L110; + } + lllen = ll_src[3]; + llfirt = 4; + goto L111; +L110: + lllen = (ll_src[5] << 15) + ll_src[4]; + llfirt = ll_src[2] + 1; +L111: + if (! (npix <= 0 || lllen <= 0)) { + goto L120; + } + ret_val = 0; + goto L100; +L120: + xe = xs + npix - 1; + skipwd = 0; + op = 1; + x1 = 1; + pv = 1; + i__1 = lllen; + for (ip = llfirt; ip <= i__1; ++ip) { + if (! skipwd) { + goto L140; + } + skipwd = 0; + goto L130; +L140: + opcode = ll_src[ip] / 4096; + data = ll_src[ip] & 4095; + sw0001 = opcode; + goto L150; +L160: + x2 = x1 + data - 1; + i1 = max(x1,xs); + i2 = min(x2,xe); + np = i2 - i1 + 1; + if (! (np > 0)) { + goto L170; + } + otop = op + np - 1; + if (! (opcode == 4)) { + goto L180; + } + i__2 = otop; + for (i__ = op; i__ <= i__2; ++i__) { + px_dst[i__] = pv; +/* L190: */ + } +/* L191: */ + goto L181; +L180: + i__2 = otop; + for (i__ = op; i__ <= i__2; ++i__) { + px_dst[i__] = 0; +/* L200: */ + } +/* L201: */ + if (! (opcode == 5 && i2 == x2)) { + goto L210; + } + px_dst[otop] = pv; +L210: +L181: + op = otop + 1; +L170: + x1 = x2 + 1; + goto L151; +L220: + pv = (ll_src[ip + 1] << 12) + data; + skipwd = 1; + goto L151; +L230: + pv += data; + goto L151; +L240: + pv -= data; + goto L151; +L250: + pv += data; + goto L91; +L260: + pv -= data; +L91: + if (! (x1 >= xs && x1 <= xe)) { + goto L270; + } + px_dst[op] = pv; + ++op; +L270: + ++x1; + goto L151; +L150: + ++sw0001; + if (sw0001 < 1 || sw0001 > 8) { + goto L151; + } + switch ((int)sw0001) { + case 1: goto L160; + case 2: goto L220; + case 3: goto L230; + case 4: goto L240; + case 5: goto L160; + case 6: goto L160; + case 7: goto L250; + case 8: goto L260; + } +L151: + if (! (x1 > xe)) { + goto L280; + } + goto L131; +L280: +L130: + ; + } +L131: + i__1 = npix; + for (i__ = op; i__ <= i__1; ++i__) { + px_dst[i__] = 0; +/* L290: */ + } +/* L291: */ + ret_val = npix; + goto L100; +L100: + return ret_val; +} /* pll2pi_ */ + Index: trunk/Ohana/src/libfits/extern/ricecomp.c =================================================================== --- trunk/Ohana/src/libfits/extern/ricecomp.c (revision 15487) +++ trunk/Ohana/src/libfits/extern/ricecomp.c (revision 15487) @@ -0,0 +1,515 @@ +/* + The following code was written by Richard White at STScI and made + available for use in CFITSIO in July 1999. These routines were + originally contained in 2 source files: rcomp.c and rdecomp.c, + and the 'include' file now called ricecomp.h was originally called buffer.h. +*/ + +/*----------------------------------------------------------*/ +/* */ +/* START OF SOURCE FILE ORIGINALLY CALLED rcomp.c */ +/* */ +/*----------------------------------------------------------*/ +/* @(#) rcomp.c 1.5 99/03/01 12:40:27 */ +/* rcomp.c Compress image line using + * (1) Difference of adjacent pixels + * (2) Rice algorithm coding + * + * Returns number of bytes written to code buffer or + * -1 on failure + */ + +#include +#include +#include +#include "ricecomp.h" /* originally included in rcomp.c file (WDP) */ +// #include "fitsio2.h" + + +static void start_outputing_bits(Buffer *buffer); +static int done_outputing_bits(Buffer *buffer); +static int output_nbits(Buffer *buffer, int bits, int n); + +/* this routine used to be called 'rcomp' (WDP) */ + +int fits_rcomp(int a[], /* input array */ + int nx, /* number of input pixels */ + unsigned char *c, /* output buffer */ + int clen, /* max length of output */ + int nblock) /* coding block size */ +{ +Buffer bufmem, *buffer = &bufmem; +int bsize, i, j, thisblock; +int lastpix, nextpix, pdiff; +int v, fs, fsmask, top, fsmax, fsbits, bbits; +int lbitbuffer, lbits_to_go; +unsigned int psum; +double pixelsum, dpsum; +unsigned int *diff; + + /* + * Original size of each pixel (bsize, bytes) and coding block + * size (nblock, pixels) + * Could make bsize a parameter to allow more efficient + * compression of short & byte images. + */ + bsize = 4; +/* nblock = 32; */ + /* + * From bsize derive: + * FSBITS = # bits required to store FS + * FSMAX = maximum value for FS + * BBITS = bits/pixel for direct coding + */ + switch (bsize) { + case 1: + fsbits = 3; + fsmax = 6; + break; + case 2: + fsbits = 4; + fsmax = 14; + break; + case 4: + fsbits = 5; + fsmax = 25; + break; + default: + fprintf (stderr, "rdecomp: bsize must be 1, 2, or 4 bytes"); + return(-1); + } + bbits = 1<start = c; + buffer->current = c; + buffer->end = c+clen; + buffer->bits_to_go = 8; + /* + * array for differences mapped to non-negative values + */ + diff = (unsigned int *) malloc(nblock*sizeof(unsigned int)); + if (diff == (unsigned int *) NULL) { + fprintf (stderr, "fits_rcomp: insufficient memory"); + return(-1); + } + /* + * Code in blocks of nblock pixels + */ + start_outputing_bits(buffer); + + /* write out first int value to the first 4 bytes of the buffer */ + if (output_nbits(buffer, a[0], 32) == EOF) { + fprintf (stderr, "rice_encode: end of buffer"); + free(diff); + return(-1); + } + + lastpix = a[0]; /* the first difference will always be zero */ + + thisblock = nblock; + for (i=0; i> 1; + for (fs = 0; psum>0; fs++) psum >>= 1; + /* + * write the codes + * fsbits ID bits used to indicate split level + */ + if (fs >= fsmax) { + /* Special high entropy case when FS >= fsmax + * Just write pixel difference values directly, no Rice coding at all. + */ + if (output_nbits(buffer, fsmax+1, fsbits) == EOF) { + fprintf (stderr, "rice_encode: end of buffer"); + free(diff); + return(-1); + } + for (j=0; jbitbuffer; + lbits_to_go = buffer->bits_to_go; + for (j=0; j> fs; + /* + * top is coded by top zeros + 1 + */ + if (lbits_to_go >= top+1) { + lbitbuffer <<= top+1; + lbitbuffer |= 1; + lbits_to_go -= top+1; + } else { + lbitbuffer <<= lbits_to_go; + if (putcbuf(lbitbuffer & 0xff,buffer) == EOF) { + fprintf (stderr, "rice_encode: end of buffer"); + free(diff); + return(-1); + } + for (top -= lbits_to_go; top>=8; top -= 8) { + if (putcbuf(0, buffer) == EOF) { + fprintf (stderr, "rice_encode: end of buffer"); + free(diff); + return(-1); + } + } + lbitbuffer = 1; + lbits_to_go = 7-top; + } + /* + * bottom FS bits are written without coding + * code is output_nbits, moved into this routine to reduce overheads + * This code potentially breaks if FS>24, so I am limiting + * FS to 24 by choice of FSMAX above. + */ + if (fs > 0) { + lbitbuffer <<= fs; + lbitbuffer |= v & fsmask; + lbits_to_go -= fs; + while (lbits_to_go <= 0) { + if (putcbuf((lbitbuffer>>(-lbits_to_go)) & 0xff,buffer)==EOF) { + fprintf (stderr, "rice_encode: end of buffer"); + free(diff); + return(-1); + } + lbits_to_go += 8; + } + } + } + buffer->bitbuffer = lbitbuffer; + buffer->bits_to_go = lbits_to_go; + } + } + done_outputing_bits(buffer); + free(diff); + /* + * return number of bytes used + */ + return(buffer->current - buffer->start); +} +/*---------------------------------------------------------------------------*/ +/* bit_output.c + * + * Bit output routines + * Procedures return zero on success, EOF on end-of-buffer + * + * Programmer: R. White Date: 20 July 1998 + */ + +/* Initialize for bit output */ + +static void start_outputing_bits(Buffer *buffer) +{ + /* + * Buffer is empty to start with + */ + buffer->bitbuffer = 0; + buffer->bits_to_go = 8; +} + +/*---------------------------------------------------------------------------*/ +/* Output N bits (N must be <= 32) */ + +static int output_nbits(Buffer *buffer, int bits, int n) +{ +/* local copies */ +int lbitbuffer; +int lbits_to_go; + + /* + * insert bits at end of bitbuffer + */ + lbitbuffer = buffer->bitbuffer; + lbits_to_go = buffer->bits_to_go; + if (lbits_to_go+n > 32) { + /* + * special case for large n: put out the top lbits_to_go bits first + * note that 0 < lbits_to_go <= 8 + */ + lbitbuffer <<= lbits_to_go; + lbitbuffer |= (bits>>(n-lbits_to_go)) & ((1<>(-lbits_to_go)) & 0xff,buffer) == EOF) + return(EOF); + lbits_to_go += 8; + } + buffer->bitbuffer = lbitbuffer; + buffer->bits_to_go = lbits_to_go; + return(0); +} + +/*---------------------------------------------------------------------------*/ +/* Flush out the last bits */ + +static int done_outputing_bits(Buffer *buffer) +{ + if(buffer->bits_to_go < 8) { + if (putcbuf(buffer->bitbuffer<bits_to_go,buffer) == EOF) + return(EOF); + } + return(0); +} +/*---------------------------------------------------------------------------*/ +/*----------------------------------------------------------*/ +/* */ +/* START OF SOURCE FILE ORIGINALLY CALLED rdecomp.c */ +/* */ +/*----------------------------------------------------------*/ + +/* @(#) rdecomp.c 1.4 99/03/01 12:38:41 */ +/* rdecomp.c Decompress image line using + * (1) Difference of adjacent pixels + * (2) Rice algorithm coding + * + * Returns 0 on success or 1 on failure + */ + +/* moved these 'includes' to the beginning of the file (WDP) +#include +#include +*/ + +/* this routine used to be called 'rdecomp' (WDP) */ + +int fits_rdecomp (unsigned char *c, /* input buffer */ + int clen, /* length of input */ + unsigned int array[], /* output array */ + int nx, /* number of output pixels */ + int nblock) /* coding block size */ +{ +int bsize, i, k, imax; +int nbits, nzero, fs; +unsigned char *cend, bytevalue; +unsigned int b, diff, lastpix; +int fsmax, fsbits, bbits; +static int *nonzero_count = (int *)NULL; + + /* + * Original size of each pixel (bsize, bytes) and coding block + * size (nblock, pixels) + * Could make bsize a parameter to allow more efficient + * compression of short & byte images. + */ + bsize = 4; +/* nblock = 32; */ + /* + * From bsize derive: + * FSBITS = # bits required to store FS + * FSMAX = maximum value for FS + * BBITS = bits/pixel for direct coding + */ + switch (bsize) { + case 1: + fsbits = 3; + fsmax = 6; + break; + case 2: + fsbits = 4; + fsmax = 14; + break; + case 4: + fsbits = 5; + fsmax = 25; + break; + default: + fprintf (stderr, "rdecomp: bsize must be 1, 2, or 4 bytes"); + return 1; + } + bbits = 1<=0; ) { + for ( ; i>=k; i--) nonzero_count[i] = nzero; + k = k/2; + nzero--; + } + } + /* + * Decode in blocks of nblock pixels + */ + + /* first 4 bytes of input buffer contain the value of the first */ + /* 4 byte integer value, without any encoding */ + + // XXX EAM : wow! this is just doing { lastpix = *(int *)c; } without a cast (encoded value is big-endian) + + lastpix = 0; + bytevalue = c[0]; + lastpix = lastpix | (bytevalue<<24); + bytevalue = c[1]; + lastpix = lastpix | (bytevalue<<16); + bytevalue = c[2]; + lastpix = lastpix | (bytevalue<<8); + bytevalue = c[3]; + lastpix = lastpix | bytevalue; + + c += 4; + cend = c + clen - 4; + + b = *c++; /* bit buffer */ + nbits = 8; /* number of bits remaining in b */ + for (i = 0; i> nbits) - 1; + b &= (1< nx) imax = nx; + if (fs<0) { + /* low-entropy case, all zero differences */ + for ( ; i= 0; k -= 8) { + b = *c++; + diff |= b<0) { + b = *c++; + diff |= b>>(-k); + b &= (1<>1; + } else { + diff = ~(diff>>1); + } + array[i] = diff+lastpix; + lastpix = array[i]; + } + } else { + /* normal case, Rice coding */ + for ( ; i>nbits); + b &= (1<>1; + } else { + diff = ~(diff>>1); + } + array[i] = diff+lastpix; + lastpix = array[i]; + } + } + if (c > cend) { + fprintf (stderr, "decompression error: hit end of compressed byte stream"); + return 1; + } + } + if (c < cend) { + fprintf (stderr, "decompression warning: unused bytes at end of compressed buffer"); + } + return 0; +} Index: trunk/Ohana/src/libfits/extern/ricecomp.h =================================================================== --- trunk/Ohana/src/libfits/extern/ricecomp.h (revision 15487) +++ trunk/Ohana/src/libfits/extern/ricecomp.h (revision 15487) @@ -0,0 +1,107 @@ +/* @(#) buffer.h 1.1 98/07/21 12:34:27 */ +/* buffer.h: structure for compression to buffer rather than to a file, including + * bit I/O buffer + * + * R. White, 19 June 1998 + */ + + +typedef unsigned char Buffer_t; + +typedef struct { + int bitbuffer; /* bit buffer */ + int bits_to_go; /* bits to go in buffer */ + Buffer_t *start; /* start of buffer */ + Buffer_t *current; /* current position in buffer */ + Buffer_t *end; /* end of buffer */ +} Buffer; + +#define buffree(mf) (free(mf->start), free(mf)) +#define bufused(mf) (mf->current - mf->start) +#define bufreset(mf) (mf->current = mf->start) + +/* + * getcbuf, putcbuf macros for character IO to buffer + * putcbuf returns EOF on end of buffer, else returns 0 + */ +#define getcbuf(mf) ((mf->current >= mf->end) ? EOF : *(mf->current)++) +#define putcbuf(c,mf) \ + ((mf->current >= mf->end) ? \ + EOF :\ + ((*(mf->current)++ = c), 0)) + +/* + * bufalloc sets up buffer of length n + */ + +/* not needed by CFITSIO + +static Buffer *bufalloc(int n) +{ +Buffer *mf; + + mf = (Buffer *) malloc(sizeof(Buffer)); + if (mf == (Buffer *)NULL) return((Buffer *)NULL); + + mf->start = (Buffer_t *) malloc(n*sizeof(Buffer_t)); + if (mf->start == (Buffer_t *)NULL) { + free(mf); + return((Buffer *)NULL); + } + mf->bits_to_go = 8; + mf->end = mf->start + n; + mf->current = mf->start; + return(mf); +} +*/ + +/* + * bufrealloc extends buffer (or truncates it) by + * reallocating memory + */ + +/* not needed by CFITSIO +static int bufrealloc(Buffer *mf, int n) +{ +int len; + + len = mf->current - mf->start; + + * silently throw away data if buffer is already longer than n * + if (len>n) len = n; + if (len<0) len = 0; + + mf->start = (Buffer_t *) realloc(mf->start, n*sizeof(Buffer_t)); + if (mf->start == (Buffer_t *)NULL) return(0); + + mf->end = mf->start + n; + mf->current = mf->start + len; + return(n); +} +*/ + +/* + * bufdump dumps contents of buffer to outfile and resets + * it to be empty. Returns number of bytes written. + * + * Note we don't write out the bit buffer -- you must call + * done_outputing_bits() first to ensure that the bit buffer + * is written out. I do it this way to allow incremental + * buffer dumps while bit IO is still going on. + */ + +/* not needed by CFITSIO + +static int bufdump(FILE *outfile, Buffer *buffer) +{ +int ndump; + + ndump = bufused(buffer); + if (fwrite(buffer->start, 1, ndump, outfile) != ndump) { + fprintf(stderr, "bufdump: error in write\n"); + exit(1); + } + bufreset(buffer); + return(ndump); +} +*/ Index: trunk/Ohana/src/libfits/header/F_convert_H.c =================================================================== --- trunk/Ohana/src/libfits/header/F_convert_H.c (revision 15487) +++ trunk/Ohana/src/libfits/header/F_convert_H.c (revision 15487) @@ -0,0 +1,52 @@ +# include +# include + +int gfits_primary_to_extended (Header *header, char *exttype, char *comment) { + + int Ns, No; + char line[81]; + + // XXX check for valid exttype, comment, output string + snprintf (line, 81, "%-8s= '%-18s' / %-s ", "XTENSION", exttype, comment); + Ns = strlen (line); + No = 80 - Ns; + strncpy (header->buffer, line, Ns); + memset (&header->buffer[Ns], ' ', No); + + return (TRUE); +} + +// don't require the current to have SIMPLE +int gfits_modify_extended (Header *header, char *exttype, char *comment) { + + int Ns, No; + char line[81]; + + // XXX check for valid exttype, comment, output string + snprintf (line, 81, "%-8s= '%-18s' / %-s ", "XTENSION", exttype, comment); + Ns = strlen (line); + No = 80 - Ns; + strncpy (header->buffer, line, Ns); + memset (&header->buffer[Ns], ' ', No); + + return (TRUE); +} + +int gfits_extended_to_primary (Header *header, int simple, char *comment) { + + int Ns, No; + char line[81]; + + // XXX check for valid exttype, comment, output string + if (simple) { + snprintf (line, 81, "%-8s= %-18s T / %-s ", "SIMPLE", " ", comment); + } else { + snprintf (line, 81, "%-8s= %-18s F / %-s ", "SIMPLE", " ", comment); + } + Ns = strlen (line); + No = 80 - Ns; + strncpy (header->buffer, line, Ns); + memset (&header->buffer[Ns], ' ', No); + + return (TRUE); +} Index: trunk/Ohana/src/libfits/header/F_copy_H.c =================================================================== --- trunk/Ohana/src/libfits/header/F_copy_H.c (revision 15240) +++ trunk/Ohana/src/libfits/header/F_copy_H.c (revision 15487) @@ -8,13 +8,18 @@ out[0].simple = in[0].simple; + out[0].unsign = in[0].unsign; + out[0].extend = in[0].extend; out[0].bitpix = in[0].bitpix; - out[0].extend = in[0].extend; - out[0].unsign = in[0].unsign; - out[0].size = in[0].size; - out[0].bscale = in[0].bscale; - out[0].bzero = in[0].bzero; + out[0].Naxes = in[0].Naxes; for (i = 0; i < FT_MAX_NAXES; i++) out[0].Naxis[i] = in[0].Naxis[i]; + + out[0].size = in[0].size; + + out[0].pcount = in[0].pcount; + out[0].gcount = in[0].gcount; + out[0].bzero = in[0].bzero; + out[0].bscale = in[0].bscale; ALLOCATE (out[0].buffer, char, out[0].size); Index: trunk/Ohana/src/libfits/header/F_create_H.c =================================================================== --- trunk/Ohana/src/libfits/header/F_create_H.c (revision 15240) +++ trunk/Ohana/src/libfits/header/F_create_H.c (revision 15487) @@ -22,11 +22,10 @@ for (i = 0; i < header[0].Naxes; i++) { - sprintf (axis, "NAXIS%d", i + 1); + snprintf (axis, 10, "NAXIS%d", i + 1); gfits_modify (header, axis, "%d", 1, header[0].Naxis[i]); } - /* PCOUNT & GCOUNT are not implemented, but they are required to be here */ - gfits_modify (header, "PCOUNT", "%d", 1, 0); - gfits_modify (header, "GCOUNT", "%d", 1, 1); + gfits_modify (header, "PCOUNT", "%d", 1, header[0].pcount); + gfits_modify (header, "GCOUNT", "%d", 1, header[0].gcount); gfits_modify (header, "BSCALE", "%lf", 1, header[0].bscale); gfits_modify (header, "BZERO", "%lf", 1, header[0].bzero); Index: trunk/Ohana/src/libfits/header/F_init_H.c =================================================================== --- trunk/Ohana/src/libfits/header/F_init_H.c (revision 15240) +++ trunk/Ohana/src/libfits/header/F_init_H.c (revision 15487) @@ -11,6 +11,8 @@ header[0].extend = FALSE; header[0].unsign = FALSE; - header[0].bscale = 1; - header[0].bzero = 0; + header[0].pcount = 0; + header[0].gcount = 1; + header[0].bscale = 1.0; + header[0].bzero = 0.0; header[0].bitpix = 8; header[0].Naxes = 0; Index: trunk/Ohana/src/libfits/header/F_read_H.c =================================================================== --- trunk/Ohana/src/libfits/header/F_read_H.c (revision 15240) +++ trunk/Ohana/src/libfits/header/F_read_H.c (revision 15487) @@ -76,4 +76,11 @@ gfits_scan (header, "NAXIS10", "%d", 1, &header[0].Naxis[9]); + if (!gfits_scan (header, "PCOUNT", "%d", 1, &header[0].pcount)) { + header[0].pcount = 0; + } + if (!gfits_scan (header, "GCOUNT", "%d", 1, &header[0].gcount)) { + header[0].gcount = 1; + } + return (TRUE); Index: trunk/Ohana/src/libfits/header/F_read_XH.c =================================================================== --- trunk/Ohana/src/libfits/header/F_read_XH.c (revision 15240) +++ trunk/Ohana/src/libfits/header/F_read_XH.c (revision 15487) @@ -45,5 +45,5 @@ } - Nmatrix = gfits_matrix_size (&theader); + Nmatrix = gfits_data_size (&theader); /* skip to next header */ Index: trunk/Ohana/src/libfits/include/gfitsio.h =================================================================== --- trunk/Ohana/src/libfits/include/gfitsio.h (revision 15240) +++ trunk/Ohana/src/libfits/include/gfitsio.h (revision 15487) @@ -1,3 +1,5 @@ /* FITS specific macros and structures */ + +# include # ifndef GFITSIO @@ -48,4 +50,6 @@ int Naxis[FT_MAX_NAXES]; int size; + int pcount; + int gcount; double bzero; double bscale; @@ -84,108 +88,118 @@ } 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 + int heap_start; // byte offset to start of HEAP + char format; // data format character (one of: XLABIJEDCM) +} VarLengthColumn; + # ifndef PROTO # define PROTO(A) A # endif -char *gfits_version PROTO(()); +char *gfits_version PROTO((void)); /******************************* Header functions *************/ -int gfits_read_header PROTO((char *, Header *)); -// int fits_copy_header PROTO((Header *, Header *)); -int gfits_copy_header PROTO((Header *, Header *)); -void gfits_free_header PROTO((Header *)); -int gfits_write_header PROTO((char *, Header *)); -int gfits_create_header PROTO((Header *)); -int gfits_delete PROTO((Header *, char *, int)); -char *gfits_header_field PROTO((Header *, char *, int)); -int gfits_read_Xheader PROTO((char *filename, Header *header, int N)); -int gfits_fread_Xheader PROTO((FILE *f, Header *header, int N)); -int gfits_load_header PROTO((FILE *f, Header *header)); - -int gfits_fread_header PROTO((FILE *, Header *)); -int gfits_fread_Theader PROTO((FILE *, Header *)); -int gfits_fread_ftable PROTO((FILE *, FTable *, char *)); -int gfits_fread_vtable_range PROTO((FILE *f, VTable *table, int start, int Nrows)); -int gfits_fread_vtable PROTO((FILE *f, VTable *table, char *extname, int Nrow, int *row)); -int gfits_fread_ftable_data PROTO((FILE *f, FTable *table)); -int gfits_fwrite_header PROTO((FILE *, Header *)); -int gfits_fwrite_Theader PROTO((FILE *, Header *)); -int gfits_fwrite_matrix PROTO((FILE *, Matrix *)); -int gfits_fwrite_table PROTO((FILE *, FTable *)); -int gfits_fwrite_vtable PROTO((FILE *, VTable *)); - -int gfits_matrix_size PROTO ((Header *)); -int gfits_init_header PROTO ((Header *)); -int gfits_get_unsign_mode PROTO (()); -int gfits_set_unsign_mode PROTO ((int mode)); -int gfits_create_table_header PROTO ((Header *, char *, char *)); -int gfits_vtable_from_ftable PROTO ((FTable *ftable, VTable *vtable, int *row, int Nrow)); -int gfits_table_to_vtable PROTO ((FTable *ftable, VTable *vtable, int start, int Nkeep)); -int gfits_create_table PROTO ((Header *, FTable *)); -int gfits_vadd_rows PROTO ((VTable *, char *, int, int)); -int gfits_add_rows PROTO ((FTable *, char *, int, int)); -char *gfits_table_print PROTO ((FTable *,...)); -int gfits_table_scale_data PROTO ((FTable *ftable)); -int gfits_table_scale_storage PROTO ((FTable *ftable)); + +char *gfits_header_field PROTO((Header *header, char *field, int N)); +char *gfits_keyword_end PROTO((char *line)); +int gfits_copy_header PROTO((Header *in, Header *out)); +int gfits_create_header PROTO((Header *header)); +int gfits_delete PROTO((Header *header, char *field, int N)); +int gfits_fread_Xheader PROTO((FILE *f, Header *header, int N)); +int gfits_fread_header PROTO((FILE *f, Header *header)); +void gfits_free_header PROTO((Header *header)); +int gfits_fwrite_header PROTO((FILE *f, Header *header)); +int gfits_get_unsign_mode PROTO((void)); +int gfits_init_header PROTO((Header *header)); +int gfits_load_header PROTO((FILE *f, Header *header)); +int gfits_modify PROTO((Header *header, char *field, char *mode, int N,...)); +int gfits_print PROTO((Header *header, char *field, char *mode, int N,...)); +int gfits_read_Xheader PROTO((char *filename, Header *header, int N)); +int gfits_read_header PROTO((char *filename, Header *header)); +int gfits_save_header PROTO((FILE *f, Header *header)); +int gfits_scan PROTO((Header *header, char *field, char *mode, int N,...)); +int gfits_set_unsign_mode PROTO((int mode)); +int gfits_stripwhite PROTO((char *string)); +int gfits_vscan PROTO((Header *header, char *field, char *mode, int N, va_list argp)); +int gfits_write_header PROTO((char *filename, Header *header)); +int gfits_data_size PROTO((Header *header)); +int gfits_extended_to_primary PROTO((Header *header, int simple, char *comment)); +int gfits_primary_to_extended PROTO((Header *header, char *exttype, char *comment)); +int gfits_modify_extended PROTO((Header *header, char *exttype, char *comment)); + /******************************* Matrix functions *************/ -int gfits_read_matrix PROTO((char *, Matrix *)); -int gfits_read_portion PROTO((char *filename, Matrix *matrix, int Nskip, int Npix)); -void gfits_free_matrix PROTO((Matrix *)); -int gfits_copy_matrix PROTO((Matrix *, Matrix *)); -int gfits_write_matrix PROTO((char *,Matrix *)); -int gfits_create_matrix PROTO((Header *, Matrix *)); -double gfits_get_matrix_value PROTO((Matrix *, int, int)); -void gfits_set_matrix_value PROTO((Matrix *, int, int, double)); -void gfits_add_matrix_value PROTO((Matrix *, int, int, double)); -int gfits_multiply_matrix PROTO((Matrix *, Matrix *, Matrix *)); -int gfits_divide_matrix PROTO((Matrix *, Matrix *, Matrix *)); -void gfits_insert_array PROTO((Matrix *, Matrix *, int, int)); -int gfits_load_matrix PROTO((FILE *f, Matrix *matrix, Header *header)); -int gfits_fread_matrix PROTO((FILE *f, Matrix *matrix, Header *header)); -int gfits_read_matrix_segment PROTO((char *filename, Matrix *matrix, char *region)); -int gfits_fread_matrix_segment PROTO((FILE *f, Matrix *matrix, Header *header, char *region)); -int gfits_convert_format PROTO((Header *header, Matrix *matrix, int outBitpix, double outScale, double outZero, int outUnsign)); + +void gfits_add_matrix_value PROTO((Matrix *matrix, int x, int y, double value)); +int gfits_convert_format PROTO((Header *header, Matrix *matrix, int outBitpix, double outScale, double outZero, int outUnsign)); +int gfits_copy_matrix PROTO((Matrix *in, Matrix *out)); +int gfits_create_matrix PROTO((Header *header, Matrix *matrix)); +int gfits_divide_matrix PROTO((Matrix *M1, Matrix *M2, Matrix *M3)); +int gfits_fread_matrix PROTO((FILE *f, Matrix *matrix, Header *header)); +int gfits_fread_matrix_segment PROTO((FILE *f, Matrix *matrix, Header *header, char *region)); +void gfits_free_matrix PROTO((Matrix *matrix)); +int gfits_fwrite_matrix PROTO((FILE *f, Matrix *matrix)); +double gfits_get_matrix_value PROTO((Matrix *matrix, int x, int y)); +void gfits_insert_array PROTO((Matrix *matrix, Matrix *array, int x, int y)); +int gfits_load_matrix PROTO((FILE *f, Matrix *matrix, Header *header)); +int gfits_multiply_matrix PROTO((Matrix *M1, Matrix *M2, Matrix *M3)); +int gfits_read_matrix PROTO((char *filename, Matrix *matrix)); +int gfits_read_matrix_segment PROTO((char *filename, Matrix *matrix, char *region)); +int gfits_read_portion PROTO((char *filename, Matrix *matrix, int Nskip, int Npix)); +void gfits_set_matrix_value PROTO((Matrix *matrix, int x, int y, double value)); +int gfits_write_matrix PROTO((char *filename, Matrix *matrix)); +int gfits_uncompress_image PROTO((Header *header, Matrix *matrix, FTable *ftable, int primary)); +int gfits_uncompress_data PROTO((char *zdata, int Nzdata, char *cmptype, char **optname, char **optvalue, int Nopt, char *outdata, int *Nout, int out_pixsize)); +int gfits_distribute_data PROTO((Matrix *matrix, int bitpix, char *data, int Ndata, int *otile, int *ztile, float zscale, float zzero)); +int gfits_byteswap_zdata PROTO((char *zdata, int Nzdata, int bitpix)); +int gfits_extension_is_compressed PROTO((Header *header)); +int gfits_tile_size PROTO((Matrix *matrix, int *otile, int *ztile)); /******************************* Table functions *************/ -int gfits_load_Theader PROTO((FILE *f, Header *Theader)); -int gfits_read_Theader PROTO((char *, Header *)); -int gfits_read_table PROTO((char *, FTable *)); -int gfits_read_table PROTO((char *, FTable *)); -int gfits_write_Theader PROTO((char *, Header *header)); -int gfits_write_table PROTO((char *, FTable *)); -int gfits_fwrite_table PROTO((FILE *f, FTable *table)); -int gfits_fwrite_vtable PROTO((FILE *f, VTable *table)); - -int gfits_bintable_format PROTO((char *, char *, int *, int *)); -int gfits_table_format PROTO((char *, char *, int *, int *)); -int gfits_save_header PROTO((FILE *, Header *header)); -int gfits_set_bintable_column PROTO((Header *header, FTable *table, char *, void *, int)); -int gfits_set_table_column PROTO((Header *header, FTable *table, char *, void *, int)); -int gfits_define_bintable_column PROTO((Header *header, char *, char *, char *, char *, double, double)); -int gfits_define_table_column PROTO((Header *header, char *, char *, char *, char *)); -int gfits_read_ftable PROTO((char *filename, FTable *table, char *extname)); - -int gfits_free_table PROTO((FTable *table)); -int gfits_free_vtable PROTO((VTable *table)); -int gfits_create_Theader PROTO((Header *header, char *type)); - -/* prototyping demanded for these four -- could be a problem on non-ANSI comp comp */ -int gfits_scan PROTO((Header *, char *, char *, int,...)); -int gfits_print PROTO((Header *, char *, char *, int,...)); -int gfits_modify PROTO((Header *, char *, char *, int,...)); -int gfits_table_column PROTO((FTable *, char *, char *,...)); -int gfits_table_column PROTO((FTable *table, char *, char *,...)); - -char *gfits_keyword_end PROTO((char *line)); -int gfits_vscan PROTO((Header *header, char *field, char *mode, int N, va_list argp)); -int gfits_scan PROTO((Header *header, char *field, char *mode, int N,...)); -int gfits_stripwhite PROTO((char *string)); -int gfits_matrix_size PROTO((Header *header)); - -int gfits_get_bintable_column_type PROTO((Header *header, char *label, char *type, int *Nval)); -int gfits_get_bintable_column PROTO((Header *header, FTable *table, char *label, void **data)); -int gfits_get_table_column_type PROTO((Header *header, char *label, char *type)); -int gfits_get_table_column PROTO((Header *header, FTable *table, char *label, void **data)); + +char *gfits_table_print PROTO((FTable *ftable,...)); +int gfits_add_rows PROTO((FTable *ftable, char *data, int Nrow, int Nbytes)); +int gfits_bintable_format PROTO((char *format, char *type, int *Nval, int *Nbytes)); +int gfits_create_table PROTO((Header *header, FTable *ftable)); +int gfits_create_table_header PROTO((Header *header, char *type, char *extname)); +int gfits_define_bintable_column PROTO((Header *header, char *format, char *label, char *comment, char *unit, double bscale, double bzero)); +int gfits_define_table_column PROTO((Header *header, char *format, char *label, char *comment, char *unit)); +int gfits_fread_ftable PROTO((FILE *f, FTable *ftable, char *extname)); +int gfits_fread_ftable_data PROTO((FILE *f, FTable *ftable)); +int gfits_fread_vtable PROTO((FILE *f, VTable *vtable, char *extname, int Nrow, int *row)); +int gfits_fread_vtable_range PROTO((FILE *f, VTable *vtable, int start, int Nrows)); +int gfits_free_table PROTO((FTable *ftable)); +int gfits_free_vtable PROTO((VTable *vtable)); +int gfits_fwrite_table PROTO((FILE *f, FTable *table)); +int gfits_fwrite_vtable PROTO((FILE *f, VTable *table)); +int gfits_get_bintable_column PROTO((Header *header, FTable *table, char *label, void **data)); +int gfits_get_bintable_column_type PROTO((Header *header, char *label, char *type, int *Nval)); +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)); +int gfits_read_ftable PROTO((char *filename, FTable *table, char *extname)); +int gfits_read_table PROTO((char *filename, FTable *ftable)); +int gfits_set_bintable_column PROTO((Header *header, FTable *table, char *label, void *data, int Nrow)); +int gfits_set_table_column PROTO((Header *header, FTable *table, char *label, void *data, int Nrow)); +int gfits_table_column PROTO((FTable *ftable, char *field, char *mode,...)); +int gfits_table_format PROTO((char *format, char *type, int *Nval, int *Nbytes)); +int gfits_table_scale_data PROTO((FTable *ftable)); +int gfits_table_scale_storage PROTO((FTable *ftable)); +int gfits_table_to_vtable PROTO((FTable *ftable, VTable *vtable, int start, int Nkeep)); +int gfits_vadd_rows PROTO((VTable *vtable, char *data, int Nrow, int Nbytes)); +int gfits_vtable_from_ftable PROTO((FTable *ftable, VTable *vtable, int *row, int Nrow)); +int gfits_write_table PROTO((char *filename, FTable *ftable)); + +int gfits_create_Theader PROTO((Header *header, char *type)); +int gfits_fread_Theader PROTO((FILE *f, Header *header)); +int gfits_fwrite_Theader PROTO((FILE *f, Header *header)); +int gfits_load_Theader PROTO((FILE *f, Header *header)); +int gfits_read_Theader PROTO((char *filename, Header *header)); +int gfits_write_Theader PROTO((char *filename, Header *header)); + +int gfits_varlength_column_define PROTO((FTable *ftable, VarLengthColumn *def, int column)); +void *gfits_varlength_column_pointer PROTO((FTable *ftable, VarLengthColumn *column, int row, int *length)); #endif /* FITSIO */ Index: trunk/Ohana/src/libfits/matrix/F_compress_M.c =================================================================== --- trunk/Ohana/src/libfits/matrix/F_compress_M.c (revision 15487) +++ trunk/Ohana/src/libfits/matrix/F_compress_M.c (revision 15487) @@ -0,0 +1,436 @@ +# include +# include +# include + +# define MOD_KEYWORD(NAME,ZNAME,TYPE,IN,OUT) { \ + gfits_scan (header, ZNAME, TYPE, 1, IN); \ + gfits_delete (header, ZNAME, 1); \ + gfits_print (header, NAME, TYPE, 1, OUT); } + +# define MOD_KEYWORD_REQUIRED(ZNAME,NAME,TYPE,IN,OUT) { \ + if (!gfits_scan (header, ZNAME, TYPE, 1, IN)) ESCAPE; \ + gfits_delete (header, ZNAME, 1); \ + gfits_print (header, NAME, TYPE, 1, OUT); } + +# define ESCAPE { \ + if (ztile != NULL) free (ztile); \ + if (optname != NULL) { \ + for (j = 0; j < Noptions; j++) { \ + free (optname[j]); \ + free (optvalue[j]); \ + } \ + free (optname); \ + free (optvalue); \ + } \ + if (out != NULL) free (out); \ + if (otile != NULL) free (otile); \ + if (ntile != NULL) free (ntile); \ + return (FALSE); } + +int gfits_uncompress_image (Header *header, Matrix *matrix, FTable *ftable, int primary) { + + int i, j, status, zimage, zcol, Nzrows, Nout, Nzdata, max_tile_size, bytes_per_pixel; + char cmptype[80]; + char zaxis[10], naxis[10], key[10], word[81], exttype[81], checksum[81], datasum[81]; + int Noptions, NOPTIONS; + VarLengthColumn zdef; + float zscale, zzero; + + int *ztile = NULL; + int *otile = NULL; + int *ntile = NULL; + char **optname = NULL; + char **optvalue = NULL; + char *out = NULL; + + // is ZIMAGE present? + status = gfits_scan (ftable->header, "ZIMAGE", "%t", 1, &zimage); + if (!status || !zimage) ESCAPE; + + // copy original header to output header + gfits_copy_header (ftable->header, header); + + // 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, "%d", &header->Naxis[i], header->Naxis[i]); + } + + // set up the tile sizes + 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 { + 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 + 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 is Noptions set to the correct value? + // XXX get ZMASKCMP + // XXX can primary be defined on basis of the header? + + if (primary) { + // XXX SIMPLE & XTENSION need to be set on the first line + if (!gfits_scan (header, "ZSIMPLE", "%t", 1, &header->simple)) ESCAPE; + gfits_delete (header, "ZSIMPLE", 1); + gfits_extended_to_primary (header, header->simple, "Image data"); + + MOD_KEYWORD ("ZEXTEND", "EXTEND", "%t", &header->extend, header->extend); + MOD_KEYWORD ("ZBLOCKED", "BLOCKED", "%t", &header->extend, header->extend); + } else { + // XXX XTENSION needs to be set on the first line + if (!gfits_scan (header, "ZTENSION", "%s", 1, exttype)) ESCAPE; + gfits_delete (header, "ZTENSION", 1); + gfits_modify_extended (header, exttype, "Image extension"); + + MOD_KEYWORD ("ZPCOUNT", "PCOUNT", "%d", &header->pcount, header->pcount); + MOD_KEYWORD ("ZGCOUNT", "GCOUNT", "%d", &header->gcount, header->gcount); + } + + MOD_KEYWORD ("ZHECKSUM", "ZHECKSUM", "%s", checksum, checksum); + MOD_KEYWORD ("ZDATASUM", "DATASUM", "%s", datasum, datasum); + + zscale = 1; + gfits_scan (header, "ZSCALE", "%lf", 1, &zscale); + + zzero = 0; + gfits_scan (header, "ZZERO", "%lf", 1, &zzero); + + // find the COMPRESSED_DATA column (format should be 1PB, 1PI, 1PJ) + for (i = 1; TRUE; i++) { + snprintf (key, 10, "TTYPE%d", i); + if (!gfits_scan (ftable->header, key, "%s", 1, word)) ESCAPE; + if (!strcmp (word, "COMPRESSED_DATA")) break; + } + zcol = i; + + if (!gfits_varlength_column_define (ftable, &zdef, zcol)) ESCAPE; + + // create the output image + gfits_create_matrix (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]; + } + + // max_tile_size is in pixels + bytes_per_pixel = abs(header[0].bitpix) / 8; + ALLOCATE (out, char, bytes_per_pixel*max_tile_size); + + // uncompress the data + Nzrows = ftable->header->Naxis[1]; + for (i = 0; i < Nzrows; i++) { + + // expected output size for this tile + Nout = gfits_tile_size (matrix, otile, ztile); + + switch (zdef.format) { + // gfits_varlength returns a pointer to a chunk of Nzdata of data elements starting at zdata. + // uncompress_data uncompresses from zdata to the temporary output buffer which must be allocated + + // XXX need to have an API to send the user data + + // XXX need to byte-swap the table column; this needs to be worked out more clearly + // in the APIs: do the table read / column extract functions swap or not?. we have + // put the byte-swapping in gifts_varlength_column_pointer, but this is fairly weak. + + case 'B': { + char *zdata; + zdata = gfits_varlength_column_pointer (ftable, &zdef, i, &Nzdata); + if (!zdata) return (FALSE); + if (!gfits_byteswap_zdata ((char *)zdata, Nzdata, 8)) return (FALSE); + if (!gfits_uncompress_data ((char *)zdata, Nzdata, cmptype, optname, optvalue, Noptions, out, &Nout, bytes_per_pixel)) return (FALSE); + if (!gfits_distribute_data (matrix, 8, out, Nout, otile, ztile, zscale, zzero)) return (FALSE); + // copy the uncompressed pixels into their correct locations + break; + } + + case 'I': { + short *zdata; + zdata = gfits_varlength_column_pointer (ftable, &zdef, i, &Nzdata); + if (!zdata) return (FALSE); + if (!gfits_byteswap_zdata ((char *)zdata, Nzdata, 16)) return (FALSE); + if (!gfits_uncompress_data ((char *)zdata, Nzdata, cmptype, optname, optvalue, Noptions, out, &Nout, bytes_per_pixel)) return (FALSE); + if (!gfits_distribute_data (matrix, 16, out, Nout, otile, ztile, zscale, zzero)) return (FALSE); + // copy the uncompressed pixels into their correct locations + break; + } + + case 'J': { + int *zdata; + if (!zdata) return (FALSE); + if (!gfits_byteswap_zdata ((char *)zdata, Nzdata, 32)) return (FALSE); + if (!gfits_uncompress_data ((char *)zdata, Nzdata, cmptype, optname, optvalue, Noptions, out, &Nout, bytes_per_pixel)) return (FALSE); + if (!gfits_distribute_data (matrix, 32, out, Nout, otile, ztile, zscale, zzero)) return (FALSE); + // copy the uncompressed pixels into their correct locations + break; + } + + default: + fprintf (stderr, "invalid size for compressed data: %c\n", zdef.format); + 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; + } + } + } + return (TRUE); +} + +int gfits_tile_size (Matrix *matrix, int *otile, int *ztile) { + + int i, Npixels, Ndimen; + + // true sizes of this tile (in pixels) + Npixels = 1; + for (i = 0; i < matrix->Naxes; i++) { + Ndimen = MIN ((matrix->Naxis[i] - otile[i]*ztile[i]), ztile[i]); + Npixels *= Ndimen; + } + + return (Npixels); +} + +// bitpix is the input data size/type +int gfits_distribute_data (Matrix *matrix, int bitpix, char *data, int Ndata, int *otile, int *ztile, float zscale, float zzero) { + + int i, j, start, offset, coord, Nline; + int *counter = NULL; + int *Ztile = NULL; + + // counter for current row in tile to copy + ALLOCATE (counter, int, matrix->Naxes); + memset (counter, 0, matrix->Naxes*sizeof(int)); + + // true sizes of this tile (in pixels) + ALLOCATE (Ztile, int, matrix->Naxes); + for (i = 0; i < matrix->Naxes; i++) { + Ztile[i] = MIN ((matrix->Naxis[i] - otile[i]*ztile[i]), ztile[i]); + } + + // number of lines in the tile (in pixels) + Nline = 1; + for (i = 1; i < matrix->Naxes; i++) { + Nline *= Ztile[i]; + } + + // double check reported size (needs the Isize) + assert (Ndata == Ztile[0]*Nline); + + // set the starting point of the tile: + // start = otile[0]*ztile[0] + otile[1]*ztile[1]*Naxis[0] + otile[2]*ztile[2]*Naxis[0]*Naxis[1] + ...; + // start = otile[0]*ztile[0] + Naxis[0]*(otile[1]*ztile[1] + Naxis[1]*(otile[2]*ztile[2] + ...)); + start = otile[matrix->Naxes-1]*ztile[matrix->Naxes-1]; + for (i = matrix->Naxes - 2; i >= 0; i--) { + coord = otile[i]*ztile[i]; + start = start*matrix->Naxis[i] + coord; + } + + // pixel offset in output array relative to tile start + offset = 0; + + // we need to set up switches for all of the possible combinations: + // this macro is used at the inner switch to run the actual loop +# define SCALE_AND_DIST(TYPE, SIZE) { \ + TYPE *Optr = (TYPE *) &matrix->buffer[SIZE*(offset + start)]; \ + for (j = 0; j < Ztile[0]; j++, Iptr++, Optr++) { \ + *Optr = *Iptr * zscale + zzero; \ + } } + + // this macro sets up the outer switch and calls above macro with all output bitpix options +# define SETUP_INSIZE(TYPE, SIZE) { \ + TYPE *Iptr = (TYPE *) &data[i*SIZE*Ztile[0]]; \ + switch (matrix->bitpix) { \ + case 8: \ + SCALE_AND_DIST (char, 1); \ + break; \ + case 16: \ + SCALE_AND_DIST (short, 2); \ + break; \ + case 32: \ + SCALE_AND_DIST (int, 4); \ + break; \ + case -32: \ + SCALE_AND_DIST (float, 4); \ + break; \ + case -64: \ + SCALE_AND_DIST (double, 8); \ + break; \ + default: \ + abort(); \ + } } + + // loop over lines in the tile + for (i = 0; i < Nline; i++) { + switch (bitpix) { + case 8: + SETUP_INSIZE (char, 1); + break; + case 16: + SETUP_INSIZE (short, 2); + break; + case 32: + SETUP_INSIZE (int, 4); + break; + default: + abort(); + } + + // update the counters, carrying to the next dimension if needed + for (j = 1; j < matrix->Naxes; j++) { + counter[j] ++; + if (counter[j] == Ztile[j]) { + counter[j] = 0; + } else { + break; + } + } + if (j == matrix->Naxes) assert (i == Nline - 1); // we should be done here... + + // Naxes = 3 + // offset = counter[1]*matrix->Naxis[0] + counter[2]*matrix->Naxis[0]*matrix->Naxis[1] + + // offset = matrix->Naxis[0]*(counter[1] + matrix->Naxis[1]*(counter[2] + matrix->Naxis[2]*...)) + + // determine the offset of the next line relative to the start position + offset = counter[matrix->Naxes - 1]; + for (j = 1; j < matrix->Naxes; j++) { + offset = offset*matrix->Naxis[j] + counter[j]; + } + } + + free (counter); + free (Ztile); + return (TRUE); +} + +int gfits_byteswap_zdata (char *zdata, int Nzdata, int bitpix) { + +# ifdef BYTE_SWAP + + int i; + char tmp; + + switch (bitpix) { + case +8: + case -8: + break; + + case +16: + case -16: + for (i = 0; i < Nzdata; i+=2) { + tmp = zdata[i]; + zdata[i] = zdata[i+1]; + zdata[i+1] = tmp; + } + break; + + case +32: + case -32: + for (i = 0; i < Nzdata; i+=4) { + tmp = zdata[i+1]; + zdata[i+1] = zdata[i+2]; + zdata[i+2] = tmp; + tmp = zdata[i]; + zdata[i] = zdata[i+3]; + zdata[i+3] = tmp; + } + break; + + case +64: + case -64: + for (i = 0; i < Nzdata; i+=8) { + tmp = zdata[i+0]; + zdata[i+0] = zdata[i+7]; + zdata[i+7] = tmp; + tmp = zdata[i+1]; + zdata[i+1] = zdata[i+6]; + zdata[i+6] = tmp; + tmp = zdata[i+2]; + zdata[i+2] = zdata[i+5]; + zdata[i+5] = tmp; + tmp = zdata[i+3]; + zdata[i+3] = zdata[i+4]; + zdata[i+4] = tmp; + } + break; + } +# endif + return (TRUE); +} + +int gfits_extension_is_compressed (Header *header) { + + int has_extname, has_zimage, zimage; + char extname[80]; + + has_extname = gfits_scan (header, "EXTNAME", "%s", 1, extname); + has_zimage = gfits_scan (header, "ZIMAGE", "%t", 1, &zimage); + + if (has_zimage && zimage) return (TRUE); + if (has_extname) { + if (!strcmp (extname, "COMPRESSED_IMAGE")) return (TRUE); + } + + return (FALSE); +} + +int gfits_compressed_is_primary (Header *header) { + + int has_ztension, has_zimage; + int ztension, zimage; + + has_zimage = gfits_scan (header, "ZIMAGE", "%t", 1, &zimage); + has_ztension = gfits_scan (header, "ZTENSION", "%t", 1, &ztension); + + if (has_zimage && zimage) return (TRUE); + + return (FALSE); +} Index: trunk/Ohana/src/libfits/matrix/F_convert_format.c =================================================================== --- trunk/Ohana/src/libfits/matrix/F_convert_format.c (revision 15240) +++ trunk/Ohana/src/libfits/matrix/F_convert_format.c (revision 15487) @@ -49,5 +49,5 @@ matrix[0].bzero = header[0].bzero = outZero; matrix[0].unsign = header[0].unsign = outUnsign; - matrix[0].size = gfits_matrix_size (header); + matrix[0].size = gfits_data_size (header); gfits_modify (header, "BITPIX", "%d", 1, outBitpix); gfits_modify (header, "BSCALE", "%lf", 1, outScale); Index: trunk/Ohana/src/libfits/matrix/F_create_M.c =================================================================== --- trunk/Ohana/src/libfits/matrix/F_create_M.c (revision 15240) +++ trunk/Ohana/src/libfits/matrix/F_create_M.c (revision 15487) @@ -18,5 +18,5 @@ matrix[0].Naxis[i] = header[0].Naxis[i]; - Nbytes = gfits_matrix_size (header); + Nbytes = gfits_data_size (header); ALLOCATE (matrix[0].buffer, char, MAX (Nbytes, 1)); bzero (matrix[0].buffer, Nbytes); Index: trunk/Ohana/src/libfits/matrix/F_get_M_value.c =================================================================== --- trunk/Ohana/src/libfits/matrix/F_get_M_value.c (revision 15240) +++ trunk/Ohana/src/libfits/matrix/F_get_M_value.c (revision 15487) @@ -3,8 +3,5 @@ /***************** fits get matrix value ***********************************/ -double gfits_get_matrix_value (matrix, x, y) -Matrix *matrix; -int x, y; -{ +double gfits_get_matrix_value (Matrix *matrix, int x, int y) { double value; Index: trunk/Ohana/src/libfits/matrix/F_insert_M.c =================================================================== --- trunk/Ohana/src/libfits/matrix/F_insert_M.c (revision 15240) +++ trunk/Ohana/src/libfits/matrix/F_insert_M.c (revision 15487) @@ -3,5 +3,5 @@ /*********************** fits insert array ***********************************/ -void gfits_insert_matrix (Matrix *matrix, Matrix *array, int x, int y) { +void gfits_insert_matrix (Matrix *matrix, Matrix *array, int x, int y) { /* there is no check here to match BITPIX, BZERO, or BSCALE */ Index: trunk/Ohana/src/libfits/matrix/F_load_M.c =================================================================== --- trunk/Ohana/src/libfits/matrix/F_load_M.c (revision 15240) +++ trunk/Ohana/src/libfits/matrix/F_load_M.c (revision 15487) @@ -25,5 +25,5 @@ matrix[0].Naxis[i] = header[0].Naxis[i]; - Nbytes = gfits_matrix_size (header); + Nbytes = gfits_data_size (header); ALLOCATE (matrix[0].buffer, char, Nbytes); matrix[0].size = Nbytes; Index: trunk/Ohana/src/libfits/matrix/F_matrix.c =================================================================== --- trunk/Ohana/src/libfits/matrix/F_matrix.c (revision 15240) +++ trunk/Ohana/src/libfits/matrix/F_matrix.c (revision 15487) @@ -2,5 +2,5 @@ # include -int gfits_matrix_size (Header *header) { +int gfits_data_size (Header *header) { int i, Nrec, size; @@ -13,4 +13,7 @@ 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) { Index: trunk/Ohana/src/libfits/matrix/F_set_M_value.c =================================================================== --- trunk/Ohana/src/libfits/matrix/F_set_M_value.c (revision 15240) +++ trunk/Ohana/src/libfits/matrix/F_set_M_value.c (revision 15487) @@ -3,9 +3,5 @@ /****************** fits set matrix value ***********************************/ -void gfits_set_matrix_value (matrix, x, y, Value) -Matrix *matrix; -int x, y; -double Value; -{ +void gfits_set_matrix_value (Matrix *matrix, int x, int y, double Value) { int pixel; Index: trunk/Ohana/src/libfits/matrix/F_uncompress_data.c =================================================================== --- trunk/Ohana/src/libfits/matrix/F_uncompress_data.c (revision 15487) +++ trunk/Ohana/src/libfits/matrix/F_uncompress_data.c (revision 15487) @@ -0,0 +1,101 @@ +# include +# include +# include + +/* functions defined in ricecomp.c */ +int fits_rcomp(int 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); + +/* 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); + +int gfits_uncompress_data (char *zdata, int Nzdata, char *cmptype, char **optname, char **optvalue, int Nopt, char *outdata, int *Nout, int out_pixsize) { + + int status; + static int Ninsum = 0; + static int Noutsum = 0; + + if (!strcasecmp(cmptype, "GZIP")) { + unsigned long tNout = *Nout; + // uncompress does not require us to know the expected number of pixel; it tells us the number + status = uncompress ((Bytef *) outdata, &tNout, (Bytef *) zdata, Nzdata); + if (status != Z_OK) { + fprintf (stderr, "error in uncompress (GZIP)\n"); + return (FALSE); + } + *Nout = tNout; + return (TRUE); + } + + if (!strcasecmp(cmptype, "RICE") || !strcasecmp(cmptype, "RICE_1")) { + 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); + } + } + } + + int Npix; + Npix = *Nout * (out_pixsize / 4.0); + + Ninsum += Nzdata; + Noutsum += *Nout; + fprintf (stderr, "%d comp bytes; %d uncomp 'pixels', totals: %d %d\n", Nzdata, Npix, Ninsum, Noutsum); + + // rice decompression from the CFITSIO source tree : we need to tell it the expected number of pixels + // is also REQUIRES 4byte output, which is fairly stupid. + if (fits_rdecomp ((unsigned char *) zdata, Nzdata, (unsigned int *) outdata, Npix, blocksize)) { + fprintf (stderr, "error in rice decompression\n"); + return (FALSE); + } + return (TRUE); + } + + if (!strcasecmp(cmptype, "PLIO")) { + int Npix; + Npix = pl_l2pi ((short *) zdata, 0, (int *) outdata, *Nout); + if (Npix != *Nout) { + fprintf (stderr, "error in plio decompression\n"); + return (FALSE); + } + return (TRUE); + } + + if (!strcasecmp(cmptype, "HCOMPRESS")) { + int Nx, Ny, scale; + status = 0; + // call hdecompress without smoothing + fits_hdecompress ((unsigned char *) zdata, FALSE, (int *) outdata, &Nx, &Ny, &scale, &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); + + if (Nx * Ny != *Nout) { + fprintf (stderr, "error in hdecompress: mismatched output size\n"); + return (FALSE); + } + return (TRUE); + } + + fprintf (stderr, "unknown compression %s\n", cmptype); + return (FALSE); +} + Index: trunk/Ohana/src/libfits/table/F_create_T.c =================================================================== --- trunk/Ohana/src/libfits/table/F_create_T.c (revision 15240) +++ trunk/Ohana/src/libfits/table/F_create_T.c (revision 15487) @@ -12,5 +12,5 @@ table[0].header = header; - Nbytes = gfits_matrix_size (header); + Nbytes = gfits_data_size (header); ALLOCATE (table[0].buffer, char, MAX (Nbytes, 1)); if (!strcmp (type, "TABLE")) { @@ -37,5 +37,5 @@ Nx = table[0].header[0].Naxis[0]; Ny = table[0].header[0].Naxis[0]; - table[0].size = gfits_matrix_size (header); + table[0].size = gfits_data_size (header); table[0].pad = table[0].size - Nx*Ny; Index: trunk/Ohana/src/libfits/table/F_create_TH.c =================================================================== --- trunk/Ohana/src/libfits/table/F_create_TH.c (revision 15240) +++ trunk/Ohana/src/libfits/table/F_create_TH.c (revision 15487) @@ -1,4 +1,5 @@ # include # include +# define NBYTES 2880 /* a basic table header (extension) is different from a primary header @@ -12,14 +13,15 @@ char axis[10]; - header[0].size = 2880; - ALLOCATE (header[0].buffer, char, 2880); + header[0].size = NBYTES; + + ALLOCATE (header[0].buffer, char, NBYTES); - for (i = 0; i < 2880; i++) + for (i = 0; i < NBYTES; i++) header[0].buffer[i] = ' '; strncpy (header[0].buffer, "END", 3); gfits_modify (header, "XTENSION", "%s", 1, type); - gfits_modify (header, "BITPIX", "%d", 1, header[0].bitpix); - gfits_modify (header, "NAXIS", "%d", 1, header[0].Naxes); + gfits_modify (header, "BITPIX", "%d", 1, header[0].bitpix); + gfits_modify (header, "NAXIS", "%d", 1, header[0].Naxes); for (i = 0; i < header[0].Naxes; i++) { @@ -28,7 +30,6 @@ } - /* PCOUNT & GCOUNT are not implemented, but they are required to be here */ - gfits_modify (header, "PCOUNT", "%d", 1, 0); - gfits_modify (header, "GCOUNT", "%d", 1, 1); + gfits_modify (header, "PCOUNT", "%d", 1, header[0].pcount); + gfits_modify (header, "GCOUNT", "%d", 1, header[0].gcount); if (!strcmp (type, "IMAGE")) { gfits_modify (header, "BSCALE", "%lf", 1, header[0].bscale); @@ -43,5 +44,4 @@ int i, valid; - char axis[10]; /* check valid table types */ @@ -51,40 +51,12 @@ if (!valid) return (FALSE); - /* create empty header buffer */ - header[0].size = 2880; - ALLOCATE (header[0].buffer, char, 2880); - - for (i = 0; i < 2880; i++) - header[0].buffer[i] = ' '; - strncpy (header[0].buffer, "END", 3); + gfits_init_header (header); - /* a table has BITPIX = 8, NAXIS = 2 */ header[0].simple = FALSE; - header[0].unsign = FALSE; - header[0].extend = FALSE; - header[0].bitpix = 8; - header[0].bscale = 1; - header[0].bzero = 0; header[0].Naxes = 2; for (i = 0; i < FT_MAX_NAXES; i++) header[0].Naxis[i] = 0; - - /* required keywords */ - gfits_modify (header, "XTENSION", "%s", 1, type); - gfits_modify (header, "BITPIX", "%d", 1, header[0].bitpix); - gfits_modify (header, "NAXIS", "%d", 1, header[0].Naxes); - - for (i = 0; i < header[0].Naxes; i++) { - sprintf (axis, "NAXIS%d", i + 1); - gfits_modify (header, axis, "%d", 1, header[0].Naxis[i], 1); - } - - /* PCOUNT & GCOUNT are not implemented, but they are required to be here */ - gfits_modify (header, "PCOUNT", "%d", 1, 0); - gfits_modify (header, "GCOUNT", "%d", 1, 1); - if (!strcmp (type, "IMAGE")) { - gfits_modify (header, "BSCALE", "%lf", 1, header[0].bscale); - gfits_modify (header, "BZERO", "%lf", 1, header[0].bzero); - } + + gfits_create_Theader (header, type); gfits_modify (header, "TFIELDS", "%d", 1, 0); Index: trunk/Ohana/src/libfits/table/F_read_T.c =================================================================== --- trunk/Ohana/src/libfits/table/F_read_T.c (revision 15240) +++ trunk/Ohana/src/libfits/table/F_read_T.c (revision 15487) @@ -41,5 +41,5 @@ /* skip to next header */ - Nbytes = gfits_matrix_size (header); + Nbytes = gfits_data_size (header); fseek (f, Nbytes, SEEK_CUR); gfits_free_header (header); @@ -54,5 +54,5 @@ /* find buffer size */ - Nbytes = gfits_matrix_size (table[0].header); + Nbytes = gfits_data_size (table[0].header); ALLOCATE (table[0].buffer, char, Nbytes); @@ -75,5 +75,5 @@ Nx = table[0].header[0].Naxis[0]; Ny = table[0].header[0].Naxis[1]; - table[0].size = gfits_matrix_size (table[0].header); + table[0].size = gfits_data_size (table[0].header); table[0].pad = table[0].size - Nx*Ny; @@ -120,5 +120,5 @@ /* find buffer size */ - Nbytes = gfits_matrix_size (header); + Nbytes = gfits_data_size (header); /* check if this is the correct extension or not */ @@ -152,64 +152,9 @@ ALLOCATE (table[0].row, int, MAX (1, Nrow)); for (i = 0; i < Nrow; i++) table[0].row[i] = row[i]; - table[0].size = gfits_matrix_size (table[0].header); + table[0].size = gfits_data_size (table[0].header); table[0].pad = table[0].size - Nx*Ny; return (TRUE); } } - -# if (0) -/*********************** fits read table ***********************************/ -/* a poor function: read only the first table, regardless of name */ -/* keep until we are sure ohana code does not use it */ -int gfits_read_table (char *filename, FTable *table) { - - FILE *f; - Header header; - int status, i, nbytes, Nbytes, matrix_size; - - if (!gfits_read_header (filename, &header)) { - fprintf (stderr, "error reading header of FITS file %s\n", filename); - return (FALSE); - } - - if (!gfits_read_Theader (filename, &table[0].header)) { - fprintf (stderr, "error reading Theader of FITS file %s\n", filename); - return (FALSE); - } - - if ((f = fopen (filename, "r")) == NULL) { - gfits_free_header (&header); - gfits_free_header (&table[0].header); - return (FALSE); - } - - matrix_size = gfits_matrix_size (&header); - - Nbytes = header.size + matrix_size + table[0].header.size; - fseek (f, Nbytes, 0); - gfits_free_header (&header); - - gfits_scan (&table[0].header, "TFIELDS", "%d", 1, &table[0].Nfields); - - /* place the field sizes in the structure? */ - table[0].bitpix = table[0].header.bitpix; - table[0].Naxes = table[0].header.Naxes; - for (i = 0; i < FT_MAX_NAXES; i++) - table[0].Naxis[i] = table[0].header.Naxis[i]; - - Nbytes = gfits_matrix_size (&table[0].header); - ALLOCATE (table[0].buffer, char, Nbytes); - nbytes = fread (table[0].buffer, sizeof(char), Nbytes, f); - if (nbytes != Nbytes) { - fprintf (stderr, "error reading in table data from FITS file %s\n", filename); - fclose (f); - return (FALSE); - } - table[0].size = Nbytes; - fclose (f); - return (TRUE); -} -# endif - int gfits_fread_header_extname (FILE *f, Header *header, char *extname) { @@ -223,5 +168,5 @@ if (!strcasecmp (extname, "PHU")) return (TRUE); - Nbytes = gfits_matrix_size (header); + Nbytes = gfits_data_size (header); fseek (f, Nbytes, SEEK_CUR); @@ -229,5 +174,5 @@ gfits_scan (header, "EXTNAME", "%s", 1, current); if (!strcmp (current, extname)) return (TRUE); - Nbytes = gfits_matrix_size (header); + Nbytes = gfits_data_size (header); fseek (f, Nbytes, SEEK_CUR); } Index: trunk/Ohana/src/libfits/table/F_read_TH.c =================================================================== --- trunk/Ohana/src/libfits/table/F_read_TH.c (revision 15240) +++ trunk/Ohana/src/libfits/table/F_read_TH.c (revision 15487) @@ -22,5 +22,5 @@ } - Nbytes = gfits_matrix_size (&header); + Nbytes = gfits_data_size (&header); fseek (f, Nbytes, SEEK_SET); gfits_free_header (&header); Index: trunk/Ohana/src/libfits/table/F_set_column.c =================================================================== --- trunk/Ohana/src/libfits/table/F_set_column.c (revision 15240) +++ trunk/Ohana/src/libfits/table/F_set_column.c (revision 15487) @@ -52,5 +52,5 @@ gfits_modify (header, "NAXIS2", "%d", 1, Ny); - nbytes = gfits_matrix_size (header); + nbytes = gfits_data_size (header); REALLOCATE (table[0].buffer, char, MAX (nbytes, 1)); bzero (table[0].buffer, nbytes); @@ -187,5 +187,5 @@ gfits_modify (header, "NAXIS2", "%d", 1, Ny); - nbytes = gfits_matrix_size (header); + nbytes = gfits_data_size (header); REALLOCATE (table[0].buffer, char, MAX (nbytes, 1)); bzero (table[0].buffer, nbytes); Index: trunk/Ohana/src/libfits/table/F_table_format.c =================================================================== --- trunk/Ohana/src/libfits/table/F_table_format.c (revision 15240) +++ trunk/Ohana/src/libfits/table/F_table_format.c (revision 15487) @@ -25,5 +25,5 @@ if (*Fchar == 'E') { *Nbytes = 4; strcpy (type, "float"); *Nval = Nv; } if (*Fchar == 'D') { *Nbytes = 8; strcpy (type, "double"); *Nval = Nv; } - if (*Fchar == 'P') { *Nbytes = 8; strcpy (type, "float"); *Nval = 2*Nv; } + if (*Fchar == 'P') { *Nbytes = 8; strcpy (type, "var"); *Nval = 2*Nv; } if (*Fchar == 'C') { *Nbytes = 8; strcpy (type, "float"); *Nval = 2*Nv; } if (*Fchar == 'M') { *Nbytes = 16; strcpy (type, "double"); *Nval = 2*Nv; } @@ -194,5 +194,5 @@ int i, j, n, Nx, Ny, Nfields; int off, Nchar, Nval, Nbytes, status; - char *line, format[64], field[16], type[16]; + char format[64], field[16], type[16]; double tzero, tscale; char *tmpChar; @@ -233,5 +233,5 @@ for (j = 0; j < Ny; j++) { for (n = 0; n < Nval; n++) { - tmpChar = (int *)&ftable[0].buffer[j*Nx + n*Nbytes + off]; + tmpChar = (char *)&ftable[0].buffer[j*Nx + n*Nbytes + off]; *tmpChar = *tmpChar * tscale + tzero; } @@ -241,5 +241,5 @@ for (j = 0; j < Ny; j++) { for (n = 0; n < Nval; n++) { - tmpShort = (int *)&ftable[0].buffer[j*Nx + n*Nbytes + off]; + tmpShort = (short *)&ftable[0].buffer[j*Nx + n*Nbytes + off]; *tmpShort = *tmpShort * tscale + tzero; } @@ -264,5 +264,5 @@ int i, j, n, Nx, Ny, Nfields; int off, Nchar, Nval, Nbytes, status; - char *line, format[64], field[16], type[16]; + char format[64], field[16], type[16]; double tzero, tscale; char *tmpChar; @@ -303,5 +303,5 @@ for (j = 0; j < Ny; j++) { for (n = 0; n < Nval; n++) { - tmpChar = (int *)&ftable[0].buffer[j*Nx + n*Nbytes + off]; + tmpChar = (char *)&ftable[0].buffer[j*Nx + n*Nbytes + off]; *tmpChar = (*tmpChar - tzero) / tscale; } @@ -311,5 +311,5 @@ for (j = 0; j < Ny; j++) { for (n = 0; n < Nval; n++) { - tmpShort = (int *)&ftable[0].buffer[j*Nx + n*Nbytes + off]; + tmpShort = (short *)&ftable[0].buffer[j*Nx + n*Nbytes + off]; *tmpShort = (*tmpShort - tzero) / tscale; } Index: trunk/Ohana/src/libfits/table/F_table_row.c =================================================================== --- trunk/Ohana/src/libfits/table/F_table_row.c (revision 15240) +++ trunk/Ohana/src/libfits/table/F_table_row.c (revision 15487) @@ -25,5 +25,5 @@ gfits_modify (header, "NAXIS2", "%d", 1, Ny); - nbytes = gfits_matrix_size (header); + nbytes = gfits_data_size (header); REALLOCATE (table[0].buffer, char, MAX (nbytes, 1)); table[0].size = nbytes; @@ -66,5 +66,5 @@ gfits_modify (header, "NAXIS2", "%d", 1, Ny); - table[0].size = gfits_matrix_size (table[0].header); + table[0].size = gfits_data_size (table[0].header); table[0].pad = table[0].size - Nx*Ny; @@ -99,5 +99,5 @@ gfits_modify (header, "NAXIS2", "%d", 1, Ny); - nbytes = gfits_matrix_size (header); + nbytes = gfits_data_size (header); REALLOCATE (table[0].buffer, char, MAX (nbytes, 1)); table[0].size = nbytes; Index: trunk/Ohana/src/libfits/table/F_table_varlength.c =================================================================== --- trunk/Ohana/src/libfits/table/F_table_varlength.c (revision 15487) +++ trunk/Ohana/src/libfits/table/F_table_varlength.c (revision 15487) @@ -0,0 +1,131 @@ +# include +# include + +int gfits_varlength_column_define (FTable *ftable, VarLengthColumn *def, int column) { + + int i, Nv, Nb; + char *p1, *p2, *p3; + char field[81]; + char format[81]; + char tmpline[81]; + + // grab the value of TFORMn for this column + snprintf (field, 80, "TFORM%d", column); + if (!gfits_scan (ftable->header, field, "%s", 1, format)) return (FALSE); + + // find and remove the max field length element + p1 = strchr (format, '('); + p2 = strchr (format, ')'); + + if (!p1 || !p2) return (FALSE); // not a valid varlength column -- missing (e_max) + if (p2 - p1 < 2) return (FALSE); // not a valid varlength column -- contains () + + def->maxlen = strtol (p1 + 1, &p3, 10); + + if (p3 != p2) return (FALSE); // not a valid varlength column -- (e_max) contains extra chars + *p1 = 0; // make the format string end here for the rest of the function + + // first char may optionally be 0, 1 + p1 = format; + if ((*p1 == '0') || (*p1 == '1')) p1 ++; + + // now p1 must be 'P'; + if (*p1 == 0) return (FALSE); + if (*p1 != 'P') return (FALSE); + + // next value is the actual varlength column format + p1 ++; + if (*p1 == 0) return (FALSE); + if (*p1 == 'P') { return (FALSE); } + + def->nbytes = 0; + if (*p1 == 'X') { def->nbytes = 1; def->format = *p1; } + if (*p1 == 'L') { def->nbytes = 1; def->format = *p1; } + if (*p1 == 'A') { def->nbytes = 1; def->format = *p1; } + if (*p1 == 'B') { def->nbytes = 1; def->format = *p1; } + if (*p1 == 'I') { def->nbytes = 2; def->format = *p1; } + if (*p1 == 'J') { def->nbytes = 4; def->format = *p1; } + if (*p1 == 'E') { def->nbytes = 4; def->format = *p1; } + if (*p1 == 'D') { def->nbytes = 8; def->format = *p1; } + if (*p1 == 'C') { def->nbytes = 8; def->format = *p1; } + if (*p1 == 'M') { def->nbytes = 16; def->format = *p1; } + if (!def->nbytes) { return (FALSE); } + + /* scan columns to find column offset */ + def->Nstart = 0; + for (i = 1; i < column; i++) { + sprintf (field, "TFORM%d", i); + gfits_scan (ftable->header, field, "%s", 1, format); + gfits_bintable_format (format, tmpline, &Nv, &Nb); + def->Nstart += Nv*Nb; + } + + if (!gfits_scan (ftable->header, "THEAP", "%d", 1, &def->heap_start)) { + def->heap_start = ftable->header->Naxis[0]*ftable->header->Naxis[1]; + } + + return TRUE; +} + +void *gfits_varlength_column_pointer (FTable *ftable, VarLengthColumn *column, int row, int *length) { + + int offset, Nx, *ptr; + + // 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 + + ptr = (int *) &ftable->buffer[row*Nx + column->Nstart]; + *length = ptr[0]; + offset = ptr[1]; + + switch (column->format) { + case 'X': + case 'L': + case 'A': + case 'B': { + char *result; + result = (char *) ftable->buffer + column->heap_start + offset; + return result; + } + case 'I': { + short *result; + result = (short *) ftable->buffer + column->heap_start + offset; + return result; + } + case 'J': { + int *result; + result = (int *) ftable->buffer + column->heap_start + offset; + return result; + } + case 'E': + case 'C': { + float *result; + result = (float *) ftable->buffer + column->heap_start + offset; + return result; + } + case 'D': + case 'M': { + double *result; + result = (double *) ftable->buffer + column->heap_start + offset; + return result; + } + case 'P': + default: + abort(); + } + abort(); +} +