/* * Compress/Uncompress a 16-bit signed int FITS image via * RICE compression. The idea is the following: * * - figure out the size of the "noisy" section of pixel values; * for example, suppose M=3 bits are considered "noisy" * * - go through the image and consider each pixel value as a * binary number * * a. take the 16-M bits that are not noisy, and convert * those bits into a decimal number. Write that * many bits, set to 0, to the output. * * b. write a bit, set to 1, to the output. * * c. write the M noisy bits to the output. * * We create an output pseudo-FITS image which has only a single * row, with a column length equal to the number of 16-bit words * in the output string (where we pad at the end to make an integral * number of 16-bit words). * * The resulting file, while it retains a FITS-like header, probably * is not standard FITS. * * Alternatively, if invoked as "unrice", then un-compress a compressed * file, turning it back into a standard 16-bit integer FITS file. * * The special "unsigned" option allows one to handle the * pseudo-FITS 16-bit unsigned integer format common in the SDSS * software. * * * * !! These programs are based upon algorithms possibly created by !! * !! Robert Rice and Pen-Shu Yeh, and definitely developed by !! * !! Nan Ellman and Robert Lupton. Nan and Robert L. have created !! * !! versions for use in dealing with Sloan Digital Sky Survey images. !! * !! All the afore-mentioned did the hard part. All I did was to !! * !! put the algorithms into practice. !! * * * * * Michael Richmond 9/07/1995 * */ #include #include #include /* needed for "strncmp" function */ #include "bits.h" /* functions in the bit-I/O source file */ #undef DEBUG #undef LINUX /* define this if you are using LINUX */ #ifndef LINUX #include #define SEEK_SET 0 /* "fseek" to move relative to beginning of file */ #define SEEK_CUR 1 /* "fseek" to move relative to current position */ #define SEEK_END 2 /* "fseek" to move relative to end of file */ #endif #define CARDLEN 80 /* # of chars in each FITS header record */ #define NCARD 36 /* # of records in a minimal header */ #define FITSMULT 2880 /* FITS file must be multiple of this in bytes */ #define SIGN_BIT 0x8000 /* the highest-order bit in a 16-bit int */ /* which we assume is the sign bit */ #define MAXDELT 128 /* if encoding would produce > this many 0 bits, */ /* then do special-case encoding instead */ #define NOISE_BITS 5 /* default value for num of noise bits */ #define MAX_ZEROES 6 /* default value for max # zeroes used in encoding */ #define FITS_EXTENSION "fts" /* new uncompressed files have this */ #define COMP_EXTENSION "ftz" /* new compressed files have this */ /* the following must be a 16-bit signed quantity */ typedef short int int16; /* and this should be a 16-bit unsigned quantity */ typedef unsigned short int uint16; /***************************************************************************** * code starts here *****************************************************************************/ /* * print the given string to the stderr, and terminate execution * with status code -1. */ static void die ( char *string /* print this on stderr */ ) { fprintf(stderr, "%s\n", string); exit(-1); } /* * compare the two string of chars, as far as the first NULL in either * string. Return * * 0 if the two are the same, up to the first NULL in the SECOND * string * 1 otherwise */ int compar ( char *s1, /* go as far as we need to in this string... */ char *s2 /* compare up to the end of this string */ ) { while ((*s1 == *s2) && (*s1 != '\0') && (*s2 != '\0')) { s1++; s2++; } if (*s2 == '\0') { return(0); } return(1); } /* * Retrieve a FITS header symbol from the header of the FILE pointer. * Return it in string form. If symbol is of type string, enclosed in * quotes, the quotes will be included as part of the string. Leading * and trailing blanks otherwise will be trimmed. * * Return 0 if we found the symbol, or 1 if not. */ int fits_get_symbol ( FILE *fits_fp, /* pointer to file */ char *symbol_name, /* name of symbol we want to find */ char *contents /* contents of the symbol's line */ ) { char line[CARDLEN + 1]; int i, j; rewind(fits_fp); while (fgets(line, CARDLEN + 1, fits_fp) != NULL) { if ((compar(line, "END") != 0) && (compar(line, symbol_name) != 0)) { continue; } if (compar(line, symbol_name) == 0) { /* skip over leading blanks */ for (i = 10; line[i] == ' '; i++) { ; } /* copy string */ for (j = 0; i < CARDLEN; j++, i++) { contents[j] = line[i]; } contents[j] = '\0'; /* remove trailing blanks */ for (j--; (contents[j] == ' ') && (j > 0); j--) { contents[j] = '\0'; } return(0); } else { return(1); } } /* we should never get here, so return an error */ return(1); } /* * Retrieve a FITS header symbol from given char-array image of * FITS header (not from a FILE pointer). * Return it in string form. If symbol is of type string, enclosed in * quotes, the quotes will be included as part of the string. Leading * and trailing blanks otherwise will be trimmed. * * Return 0 if we found the symbol, or 1 if not. */ int header_get_symbol ( char *array, /* contains image of FITS header in a file */ char *symbol_name, /* name of symbol we want to find */ char *contents /* contents of the symbol's line */ ) { char *ptr; int i, j, done; ptr = array; done = 0; while (!done) { if ((compar(ptr, "END ") != 0) && (compar(ptr, symbol_name) != 0)) { ptr += CARDLEN; continue; } if (compar(ptr, symbol_name) == 0) { /* skip over leading blanks */ for (i = 10; ptr[i] == ' '; i++) { ; } /* copy string */ for (j = 0; i < CARDLEN; j++, i++) { contents[j] = ptr[i]; } contents[j] = '\0'; /* remove trailing blanks */ for (j--; (contents[j] == ' ') && (j > 0); j--) { contents[j] = '\0'; } return(0); } else { return(1); } } /* we should never get here, so return an error */ return(1); } /* * Figure out the size of the FITS header, by looking for an "END" * card, then padding to the next-largest multiple of FITSMULT * bytes. * * Return the size of the FITS header. * */ int fits_get_header_length ( FILE *fits_fp /* pointer to file */ ) { char line[CARDLEN + 1]; long pos, length; rewind(fits_fp); while (fgets(line, CARDLEN + 1, fits_fp) != NULL) { if (compar(line, "END") == 0) { break; } } if (compar(line, "END") != 0) { die("FITS header has no END card!?"); } pos = ftell(fits_fp); length = ((pos/(long)FITSMULT) + (long)1)*(long)FITSMULT; return((int) length); } /* * Given a FILE pointer to a FITS file whose header currently * contains 'hdr_length' bytes, allocate a new char array which * is 'hdr_length' + FITSMULT bytes, and fill it with * the contents of the file's FITS header. Pad the array * with blanks. * * Return a pointer to the array, or die with error message * if we fail to allocate. */ static char * slurp_header ( FILE *fp, /* FITS file from which we'll read header */ int hdr_length /* length of this file's header, in bytes */ ) { int nbytes; char *array; nbytes = hdr_length + FITSMULT; if ((array = (char *) malloc(nbytes + 1)) == NULL) { die("slurp_header: can't allocate for header array"); } rewind(fp); if (fread(array, nbytes, 1, fp) != 1) { die("slurp_header: fread fails to read header"); } return(array); } /* * Open the given FITS file for reading. Place the number of * rows and columns into the given arguments. * Place the header length into the given argument. * Place a pointer to a (long) char array full of the header * cards into the given argument. * * Return a FILE structure for the file, or NULL if we fail to open * the file. * * If the argument "unsigned_flag" = 1, then do NOT enforce * SIMPLE = T. This flag is used to allow the convention of SDSS * unsigned 16-bit integers and SIMPLE=F. * */ FILE * fits_open ( char *fname, /* name of file to open */ int unsigned_flag, /* if 0, don't check that SIMPLE=T */ int *nrow, /* place number of row into this pointer */ int *ncol, /* place number of cols into this pointer */ int *hdr_length, /* length of FITS header, in bytes */ char **header_array /* entire FITS header, in a long string of chars */ ) { FILE *fp; char val[CARDLEN + 1]; if ((fp = fopen(fname, "r")) == NULL) { die("can't open file for reading"); } if (fits_get_symbol(fp, "SIMPLE", val) != 0) { die("fits_get_symbol fails for SIMPLE"); } if (unsigned_flag == 0) { if (compar(val, "T") != 0) { die("file is not of type SIMPLE = T"); } } if (fits_get_symbol(fp, "BITPIX", val) != 0) { die("fits_get_symbol fails for BITPIX"); } if (atoi(val) != 16) { die("file is not of type BITPIX = 16"); } if (fits_get_symbol(fp, "NAXIS", val) != 0) { die("fits_get_symbol fails for NAXIS"); } if (atoi(val) != 2) { die("file is not of type NAXIS = 2"); } if (fits_get_symbol(fp, "NAXIS1", val) != 0) { die("fits_get_symbol fails for NAXIS1"); } *ncol = atoi(val); if (fits_get_symbol(fp, "NAXIS2", val) != 0) { die("fits_get_symbol fails for NAXIS2"); } *nrow = atoi(val); *hdr_length = fits_get_header_length(fp); /* * now slurp the entire FITS header into a long character array. * we make the array long enough to accomodate one extra * set of FITSMULT bytes, beyond the present length. */ *header_array = slurp_header(fp, *hdr_length); return(fp); } /* * given an array of chars which contains a FITS header (complete with * END card and padded with blanks to a multiple of FITSMULT bytes), * write it into the given output FILE pointer. * Start writing at the very start of the output file. * */ void fits_copy_header ( char *array, /* contains image of FITS header */ FILE *out_fp /* file to which we copy FITS header */ ) { int card, nbytes, done; char *ptr; #ifdef DEBUG long pos; #endif rewind(out_fp); /* * figure out how much we have to copy to the output. What we do is * to go through the header, find the END card, and then figure * out the proper multiple of FITSMULT bytes to use. We have * to include the END card, and pad with blanks. */ ptr = array; card = 0; done = 0; while (!done) { if (compar(ptr, "END ") == 0) { done = 1; break; } ptr += CARDLEN; card++; } /* now calculate the correct number of bytes to use */ nbytes = ((card/NCARD) + 1)*NCARD*CARDLEN; /* now write out that many bytes */ if (fwrite(array, nbytes, 1, out_fp) != 1) { die("fits_copy_header: fwrite fails"); } #ifdef DEBUG pos = ftell(out_fp); printf(" pos is now %ld\n", pos); #endif } /* * pad the given FILE with null characters so that its * length is an integer multiple of FITSMULT bytes. */ void pad_file ( FILE *fp /* image file, all data + header already inside */ ) { char nulls[FITSMULT]; int i, difference; long offset, actual_length, desired_length; offset = 0; fseek(fp, offset, SEEK_END); actual_length = ftell(fp); desired_length = ((actual_length/FITSMULT) + 1)*FITSMULT; difference = desired_length - actual_length; for (i = 0; i < difference; i++) { nulls[i] = '\0'; } if (fwrite(nulls, difference, 1, fp) != 1) { die("pad_file: fwrite fails"); } } /* * Given a pointer to a position inside a big char array which * contains the image of a FITS header, and another pointer, * to a string containing a single line we want to insert into * the header, copy values from the string into the character * array. Stop copying when we hit a 'NULL' character, or * when we reach CARDLEN chars, whichever comes first. */ static void copy_fits_line ( char *insert_here, /* ptr into big array containing entire FITS header */ char *line /* ptr to line we want to insert into header */ ) { int i; for (i = 0; i < CARDLEN; i++) { if (*line != '\0') { *insert_here = *line; } else { break; } insert_here++; line++; } } /* * modify a FITS header to indicate that the data in this file is * compressed; we do so by * * 1. setting SIMPLE = F (in existing header card) * 2. setting COMPRESS= 'rice_N_K' (may require new header card) * * where N is the number of noise bits we use to compress, and * K is the maximum number of zeroes we write in an encoded number. * (so a number we DON'T encode will have K+1 zeroes). * * This, of course, makes the file no longer true FITS. */ void mark_header ( char *array, /* array full of FITS header lines */ int noise_bits, /* number of noise bits used to compress */ int max_zeroes /* max # zeroes we use to encode a number */ ) { int flag, card; char *ptr, line[CARDLEN + 2]; /* the SIMPLE= keyword ought to be the first in the file */ flag = 0; card = 0; ptr = array; while (flag == 0) { if (compar(ptr, "SIMPLE ") == 0) { sprintf(line, "SIMPLE = %20s%50s", "F", " "); copy_fits_line(ptr, line); flag = 1; break; } card++; ptr += CARDLEN; } /* * next, we add a line for "COMPRESS= 'rice_N_K'", if it isn't already * present * we place this line in the penultimate card, just before the "END" */ flag = 0; card = 0; ptr = array; while (flag == 0) { if (compar(ptr, "COMPRESS") == 0) { flag = 1; break; } if (compar(ptr, "END ") == 0) { flag = 2; break; } card++; ptr += CARDLEN; } if (flag == 0) { die("mark_header: shouldn't get here: can't find END"); } if (flag == 1) { /* * if "flag" == 1, there's already a COMPRESS line in the header; * so, let's just overwrite it. */ sprintf(line, "COMPRESS= 'rice_%02d_%03d' %56s", noise_bits, max_zeroes, " "); copy_fits_line(ptr, line); } else { /* * if we get here, "flag" = 2, meaning that no "COMPRESS" card occurs * before the current "END" card. In this case, * overwrite the current "END" card with a new COMPRESS line, * then write a new "END" card. * Note that we may safely write a new card into the array, * because we allocated extra space when we created the array. */ sprintf(line, "COMPRESS= 'rice_%02d_%03d' %56s", noise_bits, max_zeroes, " "); copy_fits_line(ptr, line); /* now write a new END card */ ptr += CARDLEN; sprintf(line, "END %76s", " "); copy_fits_line(ptr, line); } /* * remember that when we created the array, we padded it with 'blank' * characters, so no need to do that here. */ } /* * check the given char array containing the image of a FITS header * to find out if it's a pseudo-FITS file with rice-compressed data. * * Place the number of rows and columns of data into the passed args. * (although they do refer to the number of rows and cols in * UNCOMPRESSED version, not the current compressed version). * * Read the * COMPRESS = 'rice_NN_KKK' * line and parse * NN (%02d) number of noise bits used to compress * MMM (%03d) max number of zeroes used to encode * * and place them into the 'noise_bits' and 'max_zeroes' args. * * return 0 if okay * return 1 if this isn't an appropriate header. */ int check_header ( char *array, /* contains image of header from a FITS file */ int *nrow, /* number of rows in uncompressed version */ int *ncol, /* number of cols in uncompressed version */ int *noise_bits, /* number of bits used to compress */ int *max_zeroes /* max # zeroes to use when encoding a number */ ) { char val[CARDLEN + 1]; if (header_get_symbol(array, "SIMPLE", val) != 0) { fprintf(stderr, "check_header: input file has no SIMPLE keyword\n"); return(1); } if (compar(val, "F") != 0) { fprintf(stderr, "check_header: input file does not have SIMPLE = F\n"); return(1); } if (header_get_symbol(array, "COMPRESS", val) != 0) { fprintf(stderr, "check_header: input file has no COMPRESS keyword\n"); return(1); } if (compar(val, "'rice") != 0) { fprintf(stderr, "check_header: input file does not have COMPRESS= 'rice'\n"); return(1); } if (sscanf(val + 6, "%02d", noise_bits) != 1) { fprintf(stderr, "check_header: COMPRESS line says invalid number of noise bits\n"); return(1); } if (sscanf(val + 9, "%03d", max_zeroes) != 1) { fprintf(stderr, "check_header: COMPRESS line says invalid number of max_zeroes\n"); return(1); } /* get number of rows and cols */ if (header_get_symbol(array, "NAXIS1", val) != 0) { fprintf(stderr, "check_header: input file has no NAXIS1 keyword\n"); return(1); } if (sscanf(val, "%d", ncol) != 1) { fprintf(stderr, "check_header: NAXIS1 value invalid \n"); return(1); } if (header_get_symbol(array, "NAXIS2", val) != 0) { fprintf(stderr, "check_header: input file has no NAXIS2 keyword\n"); return(1); } if (sscanf(val, "%d", nrow) != 1) { fprintf(stderr, "check_header: NAXIS2 value invalid \n"); return(1); } return(0); } /* * modify a FITS header to indicate that the data in this file is * no longer compressed; we do so by * * 1. setting SIMPLE = T (unless 'unsigned_flag' == 1) * 2. setting COMPRESS= 'none' * * This, of course, makes the file FITS again. * * Place the length of the header, in bytes, into the passed arg -- * include any padding cards needed to make the header a multiple * of FITSMULT bytes. */ void unmark_header ( char *array, /* array with image of FITS header we wish */ /* to mark as "normal" */ int unsigned_flag, /* if == 1, don't set SIMPLE=T, leave as F */ int *header_length /* length of header, in bytes */ ) { int flag, card, nbytes; char *ptr, line[CARDLEN + 2]; /* the SIMPLE= keyword ought to be the first in the file */ flag = 0; ptr = array; while (flag == 0) { if (compar(ptr, "SIMPLE ") == 0) { if (unsigned_flag == 1) { sprintf(line, "SIMPLE = %20s%50s", "F", " "); } else { sprintf(line, "SIMPLE = %20s%50s", "T", " "); } copy_fits_line(ptr, line); flag = 1; break; } ptr += CARDLEN; } /* * next, we replace line * COMPRESS= 'rice_NN_KKK' * with * COMPRESS= 'none' */ flag = 0; ptr = array; card = 0; while (flag == 0) { if (compar(ptr, "COMPRESS") == 0) { flag = 1; break; } ptr += CARDLEN; card++; } /* this should never occur, of course. */ if (flag == 0) { die("unmark_header: can't find COMPRESS"); } /* over-write current COMPRESS line */ sprintf(line, "COMPRESS= 'none' %63s", " "); copy_fits_line(ptr, line); /* now continue looking down the header until we hit 'END' */ flag = 0; while (flag == 0) { if (compar(ptr, "END") == 0) { flag = 1; break; } ptr += CARDLEN; card++; } /* * calculate the header length --- must be multiple of NCARD * cards, and include the END card. */ nbytes = ((card/NCARD) + 1)*NCARD*CARDLEN; *header_length = nbytes; } /* * print the given quantity in binary (useful in debugging) */ void print_binary ( uint16 val /* quantity to be printed in binary */ ) { int i; #ifdef DEBUG printf(" val is %6u ", val); #endif for (i = 15; i >= 0; i--) { if (val & (0x1 << i)) { putchar('1'); } else { putchar('0'); } } putchar('\n'); } /* * Given a single 16-bit signed integer value, and the number * of "noisy" bits, figure out how * many bits it will take to encode it. */ int num_bits ( int16 val, /* the signed, 16-bit int value to encode */ int noisy_bits /* the number of "noise bits" at low end */ ) { uint16 uval; int total; uval = (uint16) val; #if FLIP_SIGN /* first, we create an unsigned 16-bit quantity by flipping sign bit */ if (uval & SIGN_BIT) { uval = uval & ~SIGN_BIT; } else { uval = uval | SIGN_BIT; } #endif /* next, we shift rightwards by "noisy_bits" */ uval >>= noisy_bits; #ifdef DEBUG print_binary(uval); #endif /* * the encoded version will have "uval" zeros, plus a single 1, * plus the "noisy_bits" low-order bits. * * (we haven't included special-case code for really high values yet) */ total = uval + 1 + noisy_bits; return(total); } /* * Given a single 16-bit signed integer value, and the number * of "noisy" bits, write it to the given FILE. */ int write_bits ( FILE *out_fp, /* FILE into which to write encoded data */ int16 val, /* the signed, 16-bit int value to encode */ int noisy_bits /* the number of "noise bits" at low end */ ) { uint16 uval, noise_val; uval = (uint16) val; /* this is the value of ONLY the lower-order "noisy" bits */ noise_val = uval & ~(0xFFFF << noisy_bits); /* next, we shift rightwards by "noisy_bits" */ uval >>= noisy_bits; /* * the encoded version will have * * - "uval" zeros * - a single 1, * - the "noisy_bits" low-order bits. * */ #ifdef OLD bwrite_zeroes(out_fp, uval); bwrite_ones(out_fp, 1); #else bwrite_zeroes_and_a_one(out_fp, uval); #endif #ifdef SLOW bwrite_number(out_fp, noisy_bits, noise_val); #else new_bwrite_number(out_fp, noisy_bits, noise_val); #endif return(0); } /* * Start at the given position in the file and read a set of * nrow*ncol 16-bit signed integers. Figure out the compression * factor (in bits per output "pixel") that will occur if we * use the given number of "noisy" bits. Return that compression * factor. */ float figure_compression ( FILE *in_fp, /* pointer to FITS file from which we read */ int nrow, /* number of rows of data */ int ncol, /* number of cols of data */ int header_length, /* size of FITS header, in bytes */ int noisy_bits /* treat this many low-order bits as "noise" */ ) { int nbytes, nbits, row, col, max_zeroes; int16 *buffer, *p, last_datum, delta, val; long sum; float bits_per_pix; init_buf(noisy_bits); fseek(in_fp, header_length, SEEK_SET); nbytes = 2*ncol; if ((buffer = (int16 *) malloc(nbytes)) == NULL) { die("figure_compression: can't malloc for buffer"); } /* * this is the max number of 0-bits we write for an encoded number. * * - the largest value we can send to the encoder is MAXDELT*2 * - we throw away the lowest "noisy_bits" bits from each value * * Therefore, the max number of 0-bits we can write for encoded number is * * (MAXDELT*2) / 2^(noisy_bits) = MAXDELT / 2^(noisy_bits - 1) */ max_zeroes = MAXDELT / (1 << (noisy_bits - 1)); /* and now we can read straight through the file */ last_datum = 0; sum = 0; for (row = 0; row < nrow; row++) { if (fread((char *)buffer, 2, ncol, in_fp) != ncol) { die("figure_compression: fread fails"); } p = buffer; for (col = 0; col < ncol; col++) { /* first, find difference between this datum and last one */ delta = *p - last_datum; /* * now, if abs val of difference is > MAXDELT, we write * a special set of values to the output: * * the maximum number of 0-bits an encoded value could have * one more 0-bit * a single 1-bit * 16 bits, representing 16-bit original signed pixel value * */ if ((delta > MAXDELT) || (delta < -MAXDELT)) { val = *p; nbits = max_zeroes + 1 + 1 + 16; } else { /* * but if the abs val of difference is small, then we * calculate a modified delta, and pass that to the * encoder. */ if (delta >= 0) { val = 2*delta; } else { val = (0 - delta)*2 - 1; } nbits = num_bits(val, noisy_bits); } sum += nbits; last_datum = *p++; } } /* * so the total number of compressed bits is "sum", and the total * number of pixels is "nrow*ncol"; so we can figure # bits / pixel. */ bits_per_pix = (float) sum/((float)nrow*ncol); return(bits_per_pix); } /* * Start at the given position in the file and read a set of * nrow*ncol 16-bit signed integers. Compress the data and * write an encoded version to the given output file. * * Return the average number of bits/pixel in the compressed data. */ float write_compressed_data ( FILE *in_fp, /* pointer to FITS file from which we read */ int nrow, /* number of rows of data */ int ncol, /* number of cols of data */ int header_length, /* size of FITS header, in bytes */ int noisy_bits, /* treat this many low-order bits as "noise" */ int max_zeroes, /* max. # zeroes we use to encode */ FILE *out_fp /* write output into this FILE */ ) { int nbytes, nbits, row, col; int max_delta, max_neg_delta; int16 *buffer, *p, last_datum, delta, val; long sum; float bits_per_pix; /* * set up the special bit-buffer */ init_buf(noisy_bits); fseek(in_fp, header_length, SEEK_SET); nbytes = 2*ncol; if ((buffer = (int16 *) malloc(nbytes)) == NULL) { die("write_compressed_data: can't malloc for buffer"); } /* * max_zeroes is the max number of 0-bits we write for an encoded number. * We can convert max_zeroes into the largest "delta" value that we * allow as follows: * * - we assume the lowest "noisy_bits" are all set to 1 -> 2^N - 1 * - add to that max_zeroes shifted leftwards by N bits = max_zeroes*2^N * - we then halve the sum, since this range must include both * positive and negative "delta" values * * Therefore, given max_zeroes, the largest allowed positive "delta" is * * max_delta = [ max_zeroes*2^N + (2^N - 1) ] / 2 */ max_delta = (max_zeroes*(1 << noisy_bits) + (1 << (noisy_bits - 1))) / 2; /* * we can encode the range [-(max_delta + 1) .... 0 .... max_delta] * so let's calculate the maximum negative delta here. It will save * us time in the comparison below. */ max_neg_delta = 0 - (max_delta + 1); /* and now we can read straight through the file */ last_datum = 0; sum = 0; for (row = 0; row < nrow; row++) { if (fread((char *)buffer, 2, ncol, in_fp) != ncol) { die("write_compressed_data: fread fails"); } p = buffer; for (col = 0; col < ncol; col++) { /* first, find difference between this datum and last one */ delta = *p - last_datum; /* * now, if abs val of difference is > max_delta, we write * a special set of values to the output: * * the maximum number of 0-bits an encoded value could have * one more 0-bit * a single 1-bit * 16 bits, representing 16-bit original signed pixel value * */ if ((delta > max_delta) || (delta < max_neg_delta)) { val = *p; nbits = max_zeroes + 1 + 1 + 16; #ifdef OLD bwrite_zeroes(out_fp, max_zeroes + 1); bwrite_ones(out_fp, 1); #else bwrite_zeroes_and_a_one(out_fp, max_zeroes + 1); #endif #ifdef SLOW bwrite_number(out_fp, 16, (uint16) val); #else new_bwrite_number(out_fp, 16, (uint16) val); #endif #ifdef DEBUG printf(" write (%4d, %4d) number val = %u\n", row, col, (uint16) val); #endif } else { /* * but if the abs val of difference is small, then we * calculate a modified delta, and pass that to the * encoder. */ if (delta >= 0) { val = 2*delta; } else { val = (0 - delta)*2 - 1; } #ifdef DEBUG printf(" write (%4d, %4d) delta = %d, val = %u\n", row, col, delta, val); #endif nbits = num_bits(val, noisy_bits); write_bits(out_fp, val, noisy_bits); } sum += nbits; last_datum = *p++; } } /* * flush the special bit-buffer to the file */ finish_buf(out_fp); /* * pad the file with 'null' characters to reach a multiple of FITSMULT * bytes. */ pad_file(out_fp); /* * so the total number of compressed bits is "sum", and the total * number of pixels is "nrow*ncol"; so we can figure # bits / pixel. */ bits_per_pix = (float) sum/((float)nrow*ncol); return(bits_per_pix); } /* * read a single, encoded, 16-bit unsigned int value from the given FILE. * un-encode the value and place the 16-bit SIGNED result into the * argument "val". * * return 0 if all goes well, * 1 if there is a problem */ int read_datum ( FILE *fp, /* read from this (compressed) FILE */ int max_zeroes, /* if this many 0-bits appear, no encoding */ int noisy_bits, /* number of "noisy" bits in encoded number */ int16 last_val, /* last value read from file */ uint16 *val /* place value into this argument */ ) { int nz; uint16 uval, nzeroes, noisy_val; /* * first, read the number of zeroes, plus the single 1-bit * * This is equivalent to consecutive calls to * * nz = bread_zeroes(fp); * bread_single_bit(fp); * * but it saves a function call. */ #if 0 if ((nz = bread_zeroes(fp)) == -1) { fprintf(stderr, "read_datum: can't read zeroes\n"); return(1); } nzeroes = (uint16) nz; /* next, read the single 1-bit */ if (bread_single_one(fp) != 0) { fprintf(stderr, "read_datum: can't read single 1-bit\n"); return(1); } #else if ((nz = bread_zeroes_and_a_one(fp)) == -1) { fprintf(stderr, "read_datum: can't read zeroes\n"); return(1); } nzeroes = (uint16) nz; #endif /* * now, if the number of zeroes is equal to the "maximum" number * of zeroes we can write, it means that this number has NOT * been encoded, but appears in all 16-bits following the 1-bit. * In that case, we just read the next 16 bits and set "val" * equal to their value */ if (nzeroes == max_zeroes) { if (bread_number(fp, 16, &uval) != 0) { fprintf(stderr, "read_datum: can't read 16-bit value\n"); return(1); } *val = (int16) uval; } else { /* * otherwise, the current datum IS encoded. We need to read the * next "noisy_bits" bits and then combine the number of zeroes * with the noisy_bits to build the encoded value */ if (bread_number(fp, noisy_bits, &noisy_val) != 0) { fprintf(stderr, "read_datum: can't read noisy bits\n"); return(1); } uval = (nzeroes << noisy_bits) | noisy_val; /* * now, this encoded value is the difference between the actual * pixel value and 'last_val'. If the difference is positive, * then we have * * true value = last_val + uval/2 * * whereas, if the difference is negative, we have * * true value = last_val - ((uval/2) + 1) */ if (uval % 2 == 0) { *val = last_val + uval/2; } else { *val = last_val - (uval/2 + 1); } } return(0); } /* * Start at the given position in the file and read compressed data * that should add up to nrow*ncol signed, 16-bit integer values. * write an un-encoded version to the given output file. * */ void read_compressed_data ( FILE *in_fp, /* pointer to FITS file with compressed data */ int nrow, /* number of rows of data */ int ncol, /* number of cols of data */ int header_length, /* size of FITS header of input file, in bytes */ /* should also be size of output's header */ int noisy_bits, /* treat this many low-order bits as "noise" */ int max_zeroes, /* allow this many consec. zeroes in encoded */ /* pixel value; more signals un-encoded */ FILE *out_fp /* write output into this FILE */ ) { int nbytes, row, col; int16 *buffer, *p, last_datum, val; uint16 uval; long sum; char err_msg[200]; fseek(in_fp, header_length, SEEK_SET); fseek(out_fp, header_length, SEEK_SET); /* * set up the special binary bit-buffer */ init_buf(noisy_bits); fill_buf(in_fp); nbytes = ncol*2; if ((buffer = (int16 *) malloc(nbytes)) == NULL) { die("read_compressed_data: can't malloc for buffer"); } /* and now we can read straight through the file */ last_datum = 0; sum = 0; for (row = 0; row < nrow; row++) { p = buffer; for (col = 0; col < ncol; col++) { if (read_datum(in_fp, max_zeroes + 1, noisy_bits, last_datum, &uval) != 0) { sprintf(err_msg, "read_compressed_data: error reading value for pixel (%d, %d)", row, col); die(err_msg); } val = (int16) uval; #ifdef DEBUG printf(" pixel (%4d, %4d) = %d\n", row, col, val); #endif *p++ = val; last_datum = val; } if (fwrite((char *) buffer, nbytes, 1, out_fp) != 1) { die("read_compressed_data: error writing to output file"); } } /* * now pad the output file with '\0' characters to get multiple * of FITSMULT bytes */ pad_file(out_fp); /* * free the speical bit-buffer, and the buffer we used to write * data out to the output file. */ free_buf(); free(buffer); } /* * given a file name, create a new file name which looks like * the first, but has the given extension instead of the original * extension (i.e. "foo.ftz" instead of "foo.fts"). */ char * make_extension ( char *name, /* original file name */ char *exten /* new extension for file name */ ) { int p, pos, i; char *newname; p = strlen(name); pos = p; for (p-- ; p >= 0; p--) { if (name[p] == '.') { pos = p; break; } } if ((newname = (char *) malloc(pos + strlen(exten) + 3)) == NULL) { die("make_extension: can't allocate for new name"); } strncpy(newname, name, pos); newname[pos] = '.'; for (i = 0; i < strlen(exten); i++) { newname[++pos] = exten[i]; } newname[++pos] = '\0'; #ifdef DEBUG printf("exten: old name __%s__, new name __%s__\n", name, newname); #endif return(newname); } /* * find the value of "noise_bits" that minimizes the size of the * compressed file. * * check values of noise_bits = 1 to 8. If all choices do worse * than no compression, print an error message and halt execution. */ void do_optimize ( FILE *fp, /* input file we're going to compress */ int nrow, /* number of rows in image */ int ncol, /* number of columns in image */ int header_length, /* length of FITS header, in bytes */ int *noise_bits /* the best number of noise bits */ ) { int best, bits; float min, bits_per_pix; min = 16.0; best = 0; for (bits = 1; bits < 8; bits++) { bits_per_pix = figure_compression(fp, nrow, ncol, header_length, bits); printf(" if noise_bits = %2d, then need %7.2f bits per pixel\n", bits, bits_per_pix); if (bits_per_pix < min) { best = bits; min = bits_per_pix; } } if (best == 0) { die("no compression is better than any"); } *noise_bits = best; } /* * the main program. Usage: * * rice fitsfile [optimize] [noise=N] [unsigned] * unrice compressed_fitsfile [unsigned] * * where the args are * * fitsfile FITS file to be compressed * compressed_fitsfile a file we've already rice-compressed * optimize find the best value for noise bits * by trying 1,2,...,8 * noise=N use N noise bits in compression * unsigned the original data are in pseudo-FITS * format, with SIMPLE=F, BITPIX=16 * and each pixel has unsigned 16-bit * value. Must be given both when * compressing AND un-compressing */ int main ( int argc, char *argv[] ) { char fname[100]; /* name of input file */ char out_name[100]; /* name of output file */ int i, nrow, ncol, header_length; int noise_bits; int max_zeroes; int optimize_flag, noise_flag, zeroes_flag, unsigned_flag; float bits_per_pix; int read_flag, write_flag; char *hdr_array; /* image of FITS header for output file */ FILE *in_fp; /* FILE pointer for input file */ FILE *out_fp; /* FILE pointer for output file */ read_flag = 0; write_flag = 0; optimize_flag = 0; noise_flag = 0; zeroes_flag = 0; unsigned_flag = 0; noise_bits = NOISE_BITS; /* the default value */ if (argc < 2) { fprintf(stderr, "usage: rice fitsfile [optimize] [noise=N] [unsigned]\n"); fprintf(stderr, " unrice compressed_fitsfile [unsigned]\n"); exit(-1); } /* * figure out if this program was invoked as "rice" or "unrice". * look at the last 6 characters in the zero'th argument. */ i = strlen(argv[0]); if ((i < 6) || (strncmp(argv[0] + (i - 6), "unrice", 6) != 0)) { write_flag = 1; } else { read_flag = 1; } /* parse arguments, look for optional switches */ strcpy(fname, argv[1]); for (i = 2; i < argc; i++) { if (strcmp(argv[i], "optimize") == 0) { optimize_flag = 1; continue; } if (strncmp(argv[i], "noise=", 6) == 0) { if (sscanf(argv[i] + 6, "%d", &noise_bits) != 1) { fprintf(stderr, "invalid numerical value for noise= ..%s..\n", argv[i] + 6); exit(1); } noise_flag = 1; continue; } if (strncmp(argv[i], "zeroes=", 7) == 0) { if (sscanf(argv[i] + 8, "%d", &max_zeroes) != 1) { fprintf(stderr, "invalid numerical value for zeroes= ..%s..\n", argv[i] + 6); exit(1); } zeroes_flag = 1; continue; } if (strcmp(argv[i], "unsigned") == 0) { unsigned_flag = 1; continue; } fprintf(stderr, "ignoring invalid argument ..%s..\n", argv[i]); } if ((optimize_flag == 1) && (noise_flag == 1)) { die("you must specify only one of 'optimize' and 'noise='"); } if (zeroes_flag == 0) { max_zeroes = MAX_ZEROES; } if (write_flag == 1) { in_fp = fits_open(fname, unsigned_flag, &nrow, &ncol, &header_length, &hdr_array); if (optimize_flag == 1) { do_optimize(in_fp, nrow, ncol, header_length, &noise_bits); } strcpy(out_name, make_extension(fname, COMP_EXTENSION)); out_fp = fopen(out_name, "w+"); mark_header(hdr_array, noise_bits, max_zeroes); fits_copy_header(hdr_array, out_fp); #ifdef DEBUG printf(" file %s has %d rows, %d cols, and header_length %d\n", fname, nrow, ncol, header_length); #endif } else { in_fp = fopen(fname, "r"); header_length = fits_get_header_length(in_fp); hdr_array = slurp_header(in_fp, header_length); if (check_header(hdr_array, &nrow, &ncol, &noise_bits, &max_zeroes) != 0) { die("input 'compressed' file does not have proper header"); } strcpy(out_name, make_extension(fname, FITS_EXTENSION)); out_fp = fopen(out_name, "w+"); unmark_header(hdr_array, unsigned_flag, &header_length); fits_copy_header(hdr_array, out_fp); } if (write_flag == 1) { bits_per_pix = write_compressed_data(in_fp, nrow, ncol, header_length, noise_bits, max_zeroes, out_fp); printf(" noise_bits = %2d, need %7.2f bits per pixel\n", noise_bits, bits_per_pix); } else { read_compressed_data(in_fp, nrow, ncol, header_length, noise_bits, max_zeroes, out_fp); } fclose(in_fp); fclose(out_fp); exit(0); }