Index: /branches/eam_branch_20071015/Ohana/src/libfits/Makefile =================================================================== --- /branches/eam_branch_20071015/Ohana/src/libfits/Makefile (revision 15438) +++ /branches/eam_branch_20071015/Ohana/src/libfits/Makefile (revision 15439) @@ -10,4 +10,5 @@ MATR = $(HOME)/matrix TABL = $(HOME)/table +EXT = $(HOME)/extern INC = $(HOME)/include MAN = $(HOME)/doc @@ -16,5 +17,5 @@ # programs may add their own internal requirements here FULL_CFLAGS = $(BASE_CFLAGS) -fPIC -Wall -Werror -FULL_CPPFLAGS = $(BASE_CPPFLAGS) +FULL_CPPFLAGS = $(BASE_CPPFLAGS) -I$(EXT) FULL_LDFLAGS = $(BASE_LDFLAGS) @@ -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_compress_M.$(ARCH).o +$(MATR)/F_matrix.$(ARCH).o $(MATR)/F_compress_M.$(ARCH).o \ +$(MATR)/F_uncompress_data.$(ARCH).o TABLE_OBJ = \ @@ -54,5 +56,13 @@ $(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: /branches/eam_branch_20071015/Ohana/src/libfits/extern/fits_hcompress.c =================================================================== --- /branches/eam_branch_20071015/Ohana/src/libfits/extern/fits_hcompress.c (revision 15439) +++ /branches/eam_branch_20071015/Ohana/src/libfits/extern/fits_hcompress.c (revision 15439) @@ -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: /branches/eam_branch_20071015/Ohana/src/libfits/extern/ricecomp.c =================================================================== --- /branches/eam_branch_20071015/Ohana/src/libfits/extern/ricecomp.c (revision 15439) +++ /branches/eam_branch_20071015/Ohana/src/libfits/extern/ricecomp.c (revision 15439) @@ -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: /branches/eam_branch_20071015/Ohana/src/libfits/extern/ricecomp.h =================================================================== --- /branches/eam_branch_20071015/Ohana/src/libfits/extern/ricecomp.h (revision 15439) +++ /branches/eam_branch_20071015/Ohana/src/libfits/extern/ricecomp.h (revision 15439) @@ -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: /branches/eam_branch_20071015/Ohana/src/libfits/header/F_copy_H.c =================================================================== --- /branches/eam_branch_20071015/Ohana/src/libfits/header/F_copy_H.c (revision 15438) +++ /branches/eam_branch_20071015/Ohana/src/libfits/header/F_copy_H.c (revision 15439) @@ -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: /branches/eam_branch_20071015/Ohana/src/libfits/include/gfitsio.h =================================================================== --- /branches/eam_branch_20071015/Ohana/src/libfits/include/gfitsio.h (revision 15438) +++ /branches/eam_branch_20071015/Ohana/src/libfits/include/gfitsio.h (revision 15439) @@ -89,12 +89,9 @@ typedef struct { - int length; - int offset; - int maxlen; - int nbytes; - int Nstart; - int theap; - char format; - char cmptype[81]; + 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; @@ -156,7 +153,9 @@ 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 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 *************/ Index: /branches/eam_branch_20071015/Ohana/src/libfits/matrix/F_compress_M.c =================================================================== --- /branches/eam_branch_20071015/Ohana/src/libfits/matrix/F_compress_M.c (revision 15438) +++ /branches/eam_branch_20071015/Ohana/src/libfits/matrix/F_compress_M.c (revision 15439) @@ -8,5 +8,5 @@ gfits_print (header, NAME, TYPE, 1, OUT); } -# define MOD_KEYWORD_REQUIRED(NAME,ZNAME,TYPE,IN,OUT) { \ +# define MOD_KEYWORD_REQUIRED(ZNAME,NAME,TYPE,IN,OUT) { \ if (!gfits_scan (header, ZNAME, TYPE, 1, IN)) ESCAPE; \ gfits_delete (header, ZNAME, 1); \ @@ -30,5 +30,5 @@ int gfits_uncompress_image (Header *header, Matrix *matrix, FTable *ftable, int primary) { - int i, j, status, zimage, zcol, Nzrows, Nout, Nzdata; + 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]; @@ -145,33 +145,43 @@ gfits_create_matrix (header, matrix); - // XXX is maxlen in bytes or in pixels? - ALLOCATE (out, char, zdef.maxlen); - // counters for tile numbers ALLOCATE (otile, int, matrix->Naxes); ALLOCATE (ntile, int, matrix->Naxes); + 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) { - - // may need to multiple by byte/pix - Nout = zdef.maxlen; - // gfits_varlength returns a pointer to a chunk of Nzdata of data elements starting at zdata. // uncompress_data uncompresses from zdata to the temporary output buffer which must be allocated + // XXX need to have an API to send the user data + + // 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); - gfits_byteswap_zdata (zdata, Nzdata, 8); - gfits_uncompress_data (zdata, Nzdata, zdef.cmptype, optname, optvalue, Noptions, out, &Nout); - gfits_distribute_data (matrix, 8, out, Nout, otile, ztile, zscale, zzero); // copy the uncompressed pixels into their correct locations + 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; } @@ -180,7 +190,9 @@ short *zdata; zdata = gfits_varlength_column_pointer (ftable, &zdef, i, &Nzdata); - gfits_byteswap_zdata ((char *) zdata, Nzdata, 16); - gfits_uncompress_data ((char *) zdata, Nzdata, zdef.cmptype, optname, optvalue, Noptions, out, &Nout); - gfits_distribute_data (matrix, 16, out, Nout, otile, ztile, zscale, zzero); // copy the uncompressed pixels into their correct locations + 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; } @@ -188,13 +200,15 @@ case 'J': { int *zdata; - zdata = gfits_varlength_column_pointer (ftable, &zdef, i, &Nzdata); - gfits_byteswap_zdata ((char *) zdata, Nzdata, 32); - gfits_uncompress_data ((char *) zdata, Nzdata, zdef.cmptype, optname, optvalue, Noptions, out, &Nout); - gfits_distribute_data (matrix, 32, out, Nout, otile, ztile, zscale, zzero); // copy the uncompressed pixels into their correct locations + 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: - abort(); + fprintf (stderr, "invalid size for compressed data: %c\n", zdef.format); + return (FALSE); } @@ -212,23 +226,16 @@ } -int gfits_uncompress_data (char *zdata, int Nzdata, char *cmptype, char **optname, char **optvalue, int Nopt, char *outdata, int *Nout) { - - // XXX allow lowercase value? - if (!strcasecmp(cmptype, "GZIP")) { - unsigned long tNout = *Nout; - uncompress ((Bytef *) outdata, &tNout, (Bytef *) zdata, Nzdata); - *Nout = tNout; - return (TRUE); - } - - if (!strcasecmp(cmptype, "RICE")) { +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; } - if (!strcasecmp(cmptype, "PLIO")) { - } - - if (!strcasecmp(cmptype, "HCOMPRESS")) { - } - return (TRUE); + return (Npixels); } @@ -399,2 +406,31 @@ 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: /branches/eam_branch_20071015/Ohana/src/libfits/matrix/F_uncompress_data.c =================================================================== --- /branches/eam_branch_20071015/Ohana/src/libfits/matrix/F_uncompress_data.c (revision 15439) +++ /branches/eam_branch_20071015/Ohana/src/libfits/matrix/F_uncompress_data.c (revision 15439) @@ -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: /branches/eam_branch_20071015/Ohana/src/libfits/table/F_table_varlength.c =================================================================== --- /branches/eam_branch_20071015/Ohana/src/libfits/table/F_table_varlength.c (revision 15438) +++ /branches/eam_branch_20071015/Ohana/src/libfits/table/F_table_varlength.c (revision 15439) @@ -4,5 +4,5 @@ int gfits_varlength_column_define (FTable *ftable, VarLengthColumn *def, int column) { - int i, N, Nv, Nb; + int i, Nv, Nb; char *p1, *p2, *p3; char field[81]; @@ -39,28 +39,28 @@ if (*p1 == 'P') { return (FALSE); } - def[0].nbytes = 0; - if (*p1 == 'X') { def[0].nbytes = 1; def[0].format = *p1; } - if (*p1 == 'L') { def[0].nbytes = 1; def[0].format = *p1; } - if (*p1 == 'A') { def[0].nbytes = 1; def[0].format = *p1; } - if (*p1 == 'B') { def[0].nbytes = 1; def[0].format = *p1; } - if (*p1 == 'I') { def[0].nbytes = 2; def[0].format = *p1; } - if (*p1 == 'J') { def[0].nbytes = 4; def[0].format = *p1; } - if (*p1 == 'E') { def[0].nbytes = 4; def[0].format = *p1; } - if (*p1 == 'D') { def[0].nbytes = 8; def[0].format = *p1; } - if (*p1 == 'C') { def[0].nbytes = 8; def[0].format = *p1; } - if (*p1 == 'M') { def[0].nbytes = 16; def[0].format = *p1; } - if (!def[0].nbytes) { 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[0].Nstart = 0; - for (i = 1; i < N; i++) { + 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[0].Nstart += Nv*Nb; + def->Nstart += Nv*Nb; } - if (!gfits_scan (ftable->header, "THEAP", "%d", 1, &def[0].theap)) { - def[0].theap = ftable->header->Naxis[0]*ftable->header->Naxis[1]; + if (!gfits_scan (ftable->header, "THEAP", "%d", 1, &def->heap_start)) { + def->heap_start = ftable->header->Naxis[0]*ftable->header->Naxis[1]; } @@ -70,12 +70,25 @@ void *gfits_varlength_column_pointer (FTable *ftable, VarLengthColumn *column, int row, int *length) { - int offset, Ny, *ptr; + int offset, Nx, *ptr; // find the values for the specified row // the values in the main table for this row and varlength column: - Ny = ftable->header->Naxis[1]; - ptr = (int *) &ftable->buffer[row*Ny + column->Nstart]; - offset = ptr[0]; - *length = ptr[1]; + 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) { @@ -85,15 +98,15 @@ case 'B': { char *result; - result = (char *) ftable->buffer + column->theap + offset; + result = (char *) ftable->buffer + column->heap_start + offset; return result; } case 'I': { short *result; - result = (short *) ftable->buffer + column->theap + offset; + result = (short *) ftable->buffer + column->heap_start + offset; return result; } case 'J': { int *result; - result = (int *) ftable->buffer + column->theap + offset; + result = (int *) ftable->buffer + column->heap_start + offset; return result; } @@ -101,5 +114,5 @@ case 'C': { float *result; - result = (float *) ftable->buffer + column->theap + offset; + result = (float *) ftable->buffer + column->heap_start + offset; return result; } @@ -107,5 +120,5 @@ case 'M': { double *result; - result = (double *) ftable->buffer + column->theap + offset; + result = (double *) ftable->buffer + column->heap_start + offset; return result; }