/* ------------------------------------------------------------------- Module Name: $Source: /delcano/h2/pgf/acis/huff/RCS/huff.c,v $ Purpose: Construct and test Huffman table for FITS input Assumptions: Input from 16-bit FITS file Part Number: TBD Author: Peter G. Ford References: Copyright: Massachusetts Institute of Technology 1996 $Log: huff.c,v $ * Revision 1.6 1996/07/29 17:11:14 pgf * rename 'last' to 'initval' * * Revision 1.5 1996/07/29 16:24:27 pgf * add 'initval' parameter to compressArray() and uncompressArray() * * Revision 1.4 1996/07/23 17:15:23 pgf * read and write external HUFFTAB in little-endian format * * Revision 1.3 1996/04/29 17:10:30 pgf * add tableId to HUFFTAB and -i flag to caller arguments * change previous -i flag to -r * change tableOrigin to lowLimit in HUFFTAB * add code in readTable() to check for valid lowLimit * * Revision 1.2 1996/04/27 03:49:17 pgf * add -i option and readTable() function to support it * * Revision 1.1 1996/04/17 21:33:19 pgf * Initial revision * ------------------------------------------------------------------- */ static char rcsid[] = "$Id: huff.c,v 1.6 1996/07/29 17:11:14 pgf Exp $"; /* -------------------------------------------------------------------- Function: Construct a pixel histogram from a FITS file. Then build a Huffman table and compress the file. Finally, decompress the file. ------------------------------------------------------------------- */ #include typedef struct { /* Huffman table */ unsigned tableId; /* table itentifier */ unsigned lowLimit; /* 4093 - zero difference index */ unsigned tableSize; /* number of normal entries in table */ unsigned huffTruncCode; /* tab[-3]: misc code */ unsigned huffBadBiasCode; /* tab[-2]: BAD_BIAS code */ unsigned huffBadPixelCode; /* tab[-1]: BAD_PIXEL code */ unsigned tab[1]; /* code table */ } HUFFTAB; typedef struct leaf { /* Element of Huffman tree */ struct leaf *left; /* less frequent branch */ struct leaf *right; /* more frequent branch */ unsigned freq; /* branch frequency */ int val; /* leaf value or VAL_ code */ } NODE; /* -------------------------------------------------------------------- Pixel masks and special values ------------------------------------------------------------------- */ #define PIXEL_MASK 0xfff /* pixel value */ #define BAD_PIXEL 0xfff /* ignore this pixel */ #define BAD_BIAS 0xffe /* ignore this bias */ #define HUFF_MASK 0xffffffe0 /* string area of code */ #define COUNT_MASK 0x1f /* length of code */ #define TRUNC_MAX (32-12-5) /* max length of spill */ /* -------------------------------------------------------------------- Special Huffman tree leaf values ------------------------------------------------------------------- */ #define VAL_BADPIX -1 /* bad pixel leaf val */ #define VAL_BADBIAS -2 /* bad bias leaf val */ #define VAL_MISC -3 /* out-of-range pixel */ #define VAL_NONE -4 /* leaf is a branch */ /* -------------------------------------------------------------------- Huffman table indices ------------------------------------------------------------------- */ #define HIST_MAX 8187 /* max pixel entries */ #define HIST_SIZE 8192 /* histogram size */ #define HIST_MID (HIST_MAX/2) /* center of hist[] */ #define HIST_MISC (HIST_SIZE+VAL_MISC) /* hist misc index */ #define HIST_BADBIAS (HIST_SIZE+VAL_BADBIAS) /* hist bad bias index */ #define HIST_BADPIX (HIST_SIZE+VAL_BADPIX) /* hist bad pixel index */ /* -------------------------------------------------------------------- Routines ------------------------------------------------------------------- */ unsigned short *openFits(char *, unsigned *, unsigned *); unsigned *makeHistFromFile(char *, int, unsigned); void compressFile(char *, int, char *, HUFFTAB *); void uncompressFile(char *, int, char *, HUFFTAB *); HUFFTAB *makeHuffFromFile(char *, int, unsigned, unsigned); HUFFTAB *readTable(char *); void makeTableFromTree(unsigned, unsigned, NODE *, unsigned *); int compareLeafFreqs( NODE **, NODE **); NODE *makeTreeFromTable(HUFFTAB *); NODE *addLeafToTree(NODE *, int, unsigned, unsigned); NODE *makeTreeFromHist(unsigned *, unsigned, unsigned); void insertSort(NODE **, unsigned); unsigned compressArray(unsigned short *, unsigned, unsigned *, HUFFTAB *, unsigned); unsigned uncompressArray(unsigned *, unsigned, unsigned short *, unsigned, NODE *, unsigned, unsigned); void writeTable(char *, HUFFTAB *); void swap4(unsigned *pp, unsigned count); extern char *malloc(unsigned); /*********************************************************************** Main: parse command line arguments, build trees, compress, decompress ***********************************************************************/ main(int argc, char *argv[]) { unsigned size = HIST_MAX; /* Huffman table size */ char *table = NULL; /* output table file */ int bswap = 0; /* =1 to byte-swap FITS pixels */ int init = 0; /* =1 to read table from -t */ int id = 0; /* table ID for new table */ unsigned nmisc = 2; /* number of miscellaneous counts */ HUFFTAB *huff; /* Huffman table */ while (*++argv && **argv == '-') if (argv[0][1] == 'b' && ! argv[0][2]) bswap++; else if (argv[0][1] == 'r' && ! argv[0][2]) init++; else if (argv[0][1] == 'i' && (argv[0][2] || argv[1])) id = atoi(argv[0][2] ? argv[0]+2 : *++argv); else if (argv[0][1] == 'n' && (argv[0][2] || argv[1])) size = atoi(argv[0][2] ? argv[0]+2 : *++argv); else if (argv[0][1] == 'm' && (argv[0][2] || argv[1])) nmisc = atoi(argv[0][2] ? argv[0]+2 : *++argv); else if (argv[0][1] == 't' && (argv[0][2] || argv[1])) table = argv[0][2] ? argv[0]+2 : *++argv; else fprintf(stderr, "bad -%c value\n", argv[0][1]), exit(1); if (init) { huff = readTable(table); size = huff->tableSize; } else { huff = makeHuffFromFile(argv[0], bswap, size, nmisc); huff->tableId = id; if (table) writeTable(table, huff); } compressFile(argv[0], bswap, argv[1], huff); uncompressFile(argv[0], bswap, argv[1], huff); exit(0); } /*********************************************************************** makeHuffFromFile: construct Huffman table from input file ***********************************************************************/ HUFFTAB *makeHuffFromFile(char *infile, int bswap, unsigned size, unsigned nmisc) { HUFFTAB *huff; unsigned len, maxlen = 0, jj, temp; int ii; huff = (HUFFTAB *)malloc(sizeof(HUFFTAB)+(size-1)*sizeof(unsigned)); huff->tableSize = size; huff->lowLimit = HIST_MID-size/2; makeTableFromTree(0, 0, makeTreeFromHist( makeHistFromFile(infile, bswap, nmisc), size, huff->lowLimit), huff->tab); len = huff->huffTruncCode & COUNT_MASK; if (len > TRUNC_MAX) { for (len = jj = ii = 1; ii < size; ii++) { temp = huff->tab[ii] & COUNT_MASK; if (temp <= TRUNC_MAX && temp > len) jj = ii, len = temp; } temp = huff->huffTruncCode; huff->huffTruncCode = huff->tab[jj]; huff->tab[jj] = temp; fprintf(stderr, "Warning: Huffman table rearranged %d 0x%08x <=> 0x%08x\n", jj, huff->huffTruncCode, temp); } for (ii = -3; ii < (int)size; ii++) if (maxlen < (huff->tab[ii] & COUNT_MASK)) maxlen = huff->tab[ii] & COUNT_MASK; fprintf(stderr, "Huffman %d code lengths: min %d max %d misc %d badpix %d badbias %d\n", size, huff->tab[HIST_MID-huff->lowLimit] & COUNT_MASK, maxlen, len, huff->huffBadPixelCode & COUNT_MASK, huff->huffBadBiasCode & COUNT_MASK); return huff; } /* -------------------------------------------------------------------- makeHistFromFile: construct pixel frequency histogram from 'infile' ------------------------------------------------------------------- */ unsigned *makeHistFromFile(char *infile, int bswap, unsigned nmisc) { unsigned short *inbuf; unsigned *hist, nx, ny, val, last, x, y; double sum, sumsq, npix, sqrt(double); inbuf = openFits(infile, &nx, &ny); hist = (unsigned *)malloc(HIST_SIZE*4); bzero(hist, HIST_SIZE*4); for (npix = sum = sumsq = y = 0; y < ny; y++) { if (fread(inbuf, 2, nx, stdin) != nx) perror(infile), exit(1); if (bswap) swab(inbuf, inbuf, 2*nx); last = inbuf[0] & PIXEL_MASK; if (last == BAD_PIXEL || last == BAD_BIAS) last = 0; for (x = 1; x < nx; x++) if ((val = inbuf[x] & PIXEL_MASK) == BAD_PIXEL) hist[HIST_BADPIX]++; else if (val == BAD_BIAS) hist[HIST_BADBIAS]++; else { last = val-last+HIST_MID; hist[last]++; sum += (double)last; sumsq += (double)last*(double)last; npix++; last = val; } } free(inbuf); hist[HIST_MISC] = nmisc; fprintf(stderr, "%s: input bytes %d bits %dx%dx12 mean %.2f sigma %.2f\n", infile, (3*nx*ny)/2, nx, ny, sum/npix, sqrt((sumsq-sum*sum/npix)/npix)); return hist; } /* -------------------------------------------------------------------- makeTableFromTree: build Huffman table array 'tab' from Huffman tree 'np' ------------------------------------------------------------------- */ void makeTableFromTree(unsigned code, unsigned len, NODE *np, unsigned *tab) { if (np == NULL) fprintf(stderr, "Error: null leaf pointer in makeTableFromTree\n"), exit(1); if (np->val == VAL_NONE) { makeTableFromTree(code, len+1, np->right, tab); makeTableFromTree(code+(1 << len), len+1, np->left, tab); } else { tab[np->val] = (code << (32-len)) | len; } } /* -------------------------------------------------------------------- makeTreeFromHist: build Huffman tree from pixel histogram 'hist' ------------------------------------------------------------------- */ NODE *makeTreeFromHist(unsigned *hist, unsigned size, unsigned min) { NODE *nodes[HIST_SIZE], **npp = nodes, *np; int ii; if (min < 0 || min+size > HIST_MAX) fprintf(stderr, "bad -n value: %d\n", size); /* allocate a leaf node to each pixel value */ for (ii = 0; ii < HIST_SIZE; ii++) { if (ii < min || (ii >= min+size && ii < HIST_MISC)) hist[HIST_MISC] += hist[ii]; else { np = *npp++ = (NODE *)malloc(sizeof(NODE)); np->freq = hist[ii] ? hist[ii] : 1; np->val = ii < HIST_MISC ? (ii - min) : (ii - HIST_SIZE); np->left = np->right = NULL; } } qsort(nodes, ii = npp-nodes, sizeof(NODE *), compareLeafFreqs); /* build branch nodes to connect the leaves */ for (npp = nodes, np = *npp; ii-- > 1; ) { np = (NODE *)malloc(sizeof(NODE)); np->right = *npp++; np->left = *npp; np->val = VAL_NONE; np->freq = np->right->freq + np->left->freq; *npp = np; insertSort(npp, ii); } fprintf(stderr, "Pixel frequency: max %d misc %d badpix %d badbias %d\n", hist[HIST_MID], hist[HIST_MISC], hist[HIST_BADPIX], hist[HIST_BADBIAS]); return np; } /* -------------------------------------------------------------------- compareLeafFreqs: called from qsort() to compare two leaf frequencies ------------------------------------------------------------------- */ int compareLeafFreqs( NODE **a, NODE **b) { return a[0]->freq - b[0]->freq; } /* -------------------------------------------------------------------- insertSort: insertion sort of node npp[0] into npp[] array ------------------------------------------------------------------- */ void insertSort(NODE **npp, unsigned ii) { NODE *np; for (np = npp[0]; ii-- > 1 && np->freq > npp[1]->freq; npp++) npp[0] = npp[1]; npp[0] = np; } /*********************************************************************** compressFile: compress 'infile' to 'outfile' using Huffman table 'huff' ***********************************************************************/ void compressFile(char *infile, int bswap, char *outfile, HUFFTAB *huff) { FILE *fp; unsigned nx, ny, y, *outbuf, nw; unsigned short *inbuf, outlen; inbuf = openFits(infile, &nx, &ny); outbuf = (unsigned *)malloc((nx+1)*4); if ((fp = fopen(outfile, "w")) == NULL) perror(outfile), exit(1); fwrite(&nx, 1, 4, fp); fwrite(&ny, 1, 4, fp); for (y = nw = 0; y < ny; y++) { fread(inbuf, nx, 2, stdin); if (bswap) swab(inbuf, inbuf, 2*nx); outlen = compressArray(inbuf, nx, outbuf, huff, 0); fwrite(&outlen, 1, 2, fp); fwrite(outbuf, outlen, 4, fp); nw += outlen; } fclose(fp); free(inbuf); free(outbuf); fprintf(stderr, "%s: compressed to %d bytes (%.2f%%)\n", outfile, 4*nw+8, (800.0*nw)/(3*nx*ny)); } /* -------------------------------------------------------------------- compressArray: use Huffman table 'huff' to compress 'in' array to 'out' ------------------------------------------------------------------- */ unsigned compressArray( unsigned short *in, /* input pixel array */ unsigned inlen, /* number of input pixels */ unsigned *out, /* output array */ HUFFTAB *huff, /* Huffman table */ unsigned initval /* previous pixel value */ ) { unsigned *outorg = out; /* saved output origin */ unsigned val; /* input pixel value */ unsigned code; /* coded pixel value */ unsigned bitlen = 12; /* length of 'code' in bits */ unsigned acc = 0; /* output register */ unsigned bitout = 0; /* length of 'acc' in bits */ unsigned trunc; /* truncated pixel code */ int index; /* Huffman table index */ /* construct code for truncated pixels */ trunc = (12 + (huff->huffTruncCode & COUNT_MASK)) | ((huff->huffTruncCode & HUFF_MASK) >> 12); /* compress the array */ while (inlen--) { if (huff) { /* load the next pixel */ switch (val = *in++ & PIXEL_MASK) { case BAD_PIXEL: code = huff->huffBadPixelCode; break; case BAD_BIAS: code = huff->huffBadBiasCode; break; default: index = (val+(HIST_MID-huff->lowLimit))-initval; initval = val; if (index < 0 || index >= huff->tableSize) code = trunc | (val << 20); else code = huff->tab[index]; break; } bitlen = code & COUNT_MASK; code >>= 32 - bitlen; } else code = *in++; /* append 'code' to 'acc' */ acc |= code << bitout; if ((bitout += bitlen ) >= 32) { bitout -= 32; *out++ = acc; acc = code >> (bitlen - bitout); } } /* anything left to save? */ if (bitout > 0) *out++ = acc; /* return output word count */ return out - outorg; } /*********************************************************************** uncompressFile: decompress 'tmpfile' using Huffman tree 'root' and compare the result to the original 'infile' ***********************************************************************/ void uncompressFile(char *infile, int bswap, char *tmpfile, HUFFTAB *huff) { FILE *fp; NODE *root; unsigned nx, ny, x, y, ii, maxlen, *tmpbuf = NULL; unsigned short *inbuf, *cmpbuf, tmplen; root = makeTreeFromTable(huff); inbuf = openFits(infile, &nx, &ny); if ((fp = fopen(tmpfile, "r")) == NULL) perror(tmpfile), exit(1); if (fread(&ii, 4, 1, fp) != 1 || ii != nx) fprintf(stderr, "%s: bad x-dim: %d (not %d)\n", tmpfile, ii, nx); if (fread(&ii, 4, 1, fp) != 1 || ii != ny) fprintf(stderr, "%s: bad y-dim: %d (not %d)\n", tmpfile, ii, ny); cmpbuf = (unsigned short *)malloc(nx*2); for (y = maxlen = 0; y < ny; y++) { if (fread(inbuf, 2, nx, stdin) != nx) perror(infile), exit(1); if (bswap) swab(inbuf, inbuf, nx*2); if (fread(&tmplen, 2, 1, fp) != 1) perror(tmpfile), exit(1); if (tmplen > maxlen) { if (tmpbuf) free(tmpbuf); tmpbuf = (unsigned *)malloc((maxlen = tmplen)*4); } if (fread(tmpbuf, 4, tmplen, fp) != tmplen) perror(tmpfile), exit(1); if (!uncompressArray(tmpbuf, tmplen, cmpbuf, nx, root, HIST_MID-huff->lowLimit, 0)) fprintf(stderr, "%s: decompression fails at line %d\n", tmpfile, y), exit(1); for (x = 0; x < nx; x++) if ((inbuf[x] & PIXEL_MASK) != cmpbuf[x]) fprintf(stderr, "%s: line %4d col %4d orig %5d unpk %5d\n", tmpfile, y, x, inbuf[x], cmpbuf[x]); } fclose(fp); free(inbuf); free(tmpbuf); free(cmpbuf); fprintf(stderr, "%s: uncompressFile was successful\n", tmpfile); } /* -------------------------------------------------------------------- makeTreeFromTable: construct a Huffman tree from Huffman table 'huff' ------------------------------------------------------------------- */ NODE *makeTreeFromTable(HUFFTAB *huff) { NODE *root = NULL; unsigned tab[HIST_SIZE]; int ii; for (ii = -3; ii < (int)huff->tableSize; ii++) root = addLeafToTree(root, ii, huff->tab[ii] & HUFF_MASK, huff->tab[ii] & COUNT_MASK); bzero(tab, sizeof(tab)); makeTableFromTree(0, 0, root, tab+3); if (bcmp(tab, &huff->huffTruncCode, 4*(huff->tableSize+3))) fprintf(stderr, "failed to construct Huffman tree\n"); return root; } /* -------------------------------------------------------------------- addLeafToTree: add a new leaf to a Huffman tree ------------------------------------------------------------------- */ NODE *addLeafToTree(NODE *np, int val, unsigned code, unsigned len) { if (np == NULL) { np = (NODE *)malloc(sizeof(NODE)); np->left = np->right = NULL; np->freq = 0; np->val = VAL_NONE; } if (len == 0) np->val = val; else if (code & (1 << (32 - len))) np->left = addLeafToTree(np->left, val, code, len-1); else np->right = addLeafToTree(np->right, val, code, len-1); return np; } /* -------------------------------------------------------------------- uncompressArray: use Huffman tree 'root' to unpack array 'in' to 'out' ------------------------------------------------------------------- */ unsigned uncompressArray( unsigned *in, /* input packed array */ unsigned inlen, /* number of input elements */ unsigned short *out, /* output array */ unsigned outlen, /* maximum length of output */ NODE *root, /* Huffman tree */ unsigned offset, /* Huffman table zero difference index */ unsigned initval /* previous unpacked value */ ) { unsigned code; /* unpacking buffer */ unsigned bitlen; /* length of unpacking buffer */ unsigned temp; /* scratch */ NODE *np = root; /* Huffman tree pointer */ while (inlen-- > 0) { code = *in++; bitlen = 32; while (bitlen--) { np = (code & 1) ? np->left : np->right; code >>= 1; switch (np->val) { case VAL_NONE: continue; case VAL_BADPIX: *out++ = BAD_PIXEL; break; case VAL_BADBIAS: *out++ = BAD_BIAS; break; case VAL_MISC: if (bitlen >= 12) { initval = *out++ = code & 0xfff; code >>= 12; bitlen -= 12; } else { temp = code; if (inlen-- == 0) return 0; code = *in++; initval = *out++ = (temp | (code << bitlen)) & 0xfff; bitlen += 20; code >>= 32-bitlen; } break; default: initval = *out++ = np->val+initval-offset; break; } if (--outlen == 0) return inlen == 0; np = root; } } return 0; } /*********************************************************************** openFits: open a FITS file using stdin, read its header, return nx,ny ***********************************************************************/ unsigned short *openFits(char *infile, unsigned *nx, unsigned *ny) { char hdr[2880]; int ii, neof, bits = 0; if (freopen(infile, "r", stdin) == NULL) perror(infile), exit(1); *nx = *ny = 0; do { if (fread(hdr, sizeof(hdr), 1, stdin) != 1) perror(infile), exit(1); for (ii = 0; ii < 2880 && (neof = strncmp(hdr+ii, "END ", 4)); ii += 80) { sscanf(hdr+ii, "NAXIS1 = %d", nx); sscanf(hdr+ii, "NAXIS2 = %d", ny); sscanf(hdr+ii, "BITPIX = %d", &bits); } } while(neof); if (bits != 16 || *nx <= 0 || *ny <= 0) fprintf(stderr, "%s: bad FITS header\n", infile), exit(1); return (unsigned short *)malloc(*nx*2); } /*********************************************************************** readTable: read Huffman table from file 'file' ***********************************************************************/ HUFFTAB *readTable(char *table) { FILE *fp; HUFFTAB *huff; unsigned size; if ((fp = fopen(table, "r")) == NULL || fseek(fp, 0L, 2)) perror(table), exit(1); huff = (HUFFTAB *)malloc(size = ftell(fp)); if (fseek(fp, 0L, 0) || fread(huff, size, 1, fp) != 1) perror(table), exit(1); (void)fclose(fp); swap4((unsigned *)huff, size/4); fprintf(stderr, "%s: id %d size %d low %d\n", table, huff->tableId, huff->tableSize, huff->lowLimit); if (huff->lowLimit+huff->tableSize >= HIST_MAX) fprintf(stderr, "%s: bad lowLimit value\n", table), exit(1); return huff; } /*********************************************************************** writeTable: write Huffman table 'huff' to file 'outfile' ***********************************************************************/ void writeTable(char *table, HUFFTAB *huff) { FILE *fp; unsigned size = sizeof(HUFFTAB)+(huff->tableSize-1)*sizeof(unsigned); if ((fp = fopen(table, "w")) == NULL) perror(table), exit(1); swap4((unsigned *)huff, size/4); if (fwrite(huff, size, 1, fp) != 1) perror(table), exit(1); (void)fclose(fp); swap4((unsigned *)huff, size/4); } /*********************************************************************** swap4: perform byte-reversal of 32-bit array ***********************************************************************/ void swap4(unsigned *pp, unsigned count) { #ifndef mips #ifndef vax while (count-- > 0) { int ii = *pp; *pp++ = ((ii & 0xff) << 24) | ((ii & 0xff00) << 8) | ((ii >> 8) & 0xff00) | ((ii >> 24) & 0xff); } #endif vax #endif mips } /* -------------------------------------------------------------------- End ------------------------------------------------------------------- */