aboutsummaryrefslogtreecommitdiff
path: root/vendor/cfitsio/fits_hcompress.c
diff options
context:
space:
mode:
Diffstat (limited to 'vendor/cfitsio/fits_hcompress.c')
-rw-r--r--vendor/cfitsio/fits_hcompress.c1858
1 files changed, 1858 insertions, 0 deletions
diff --git a/vendor/cfitsio/fits_hcompress.c b/vendor/cfitsio/fits_hcompress.c
new file mode 100644
index 00000000..96a6b126
--- /dev/null
+++ b/vendor/cfitsio/fits_hcompress.c
@@ -0,0 +1,1858 @@
+/* #########################################################################
+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 <stdio.h>
+#include <string.h>
+#include <math.h>
+#include <stdlib.h>
+#include "fitsio2.h"
+
+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 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 qwrite(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(void);
+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) */
+static void output_nybble(char *outfile, int bits);
+static void output_nnybble(char *outfile, int n, unsigned char array[]);
+
+#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 */
+
+ FFLOCK;
+ noutmax = *nbytes; /* input value is the allocated size of the array */
+ *nbytes = 0; /* reset */
+
+ stat = encode(output, nbytes, a, nx, ny, scale);
+ FFUNLOCK;
+
+ *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 */
+
+ FFLOCK;
+ noutmax = *nbytes; /* input value is the allocated size of the array */
+ *nbytes = 0; /* reset */
+
+ stat = encode64(output, nbytes, a, nx, ny, scale);
+ FFUNLOCK;
+
+ *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<<log2n) ) {
+ log2n += 1;
+ }
+ /*
+ * get temporary storage for shuffling elements
+ */
+ tmp = (int *) malloc(((nmax+1)/2)*sizeof(int));
+ if(tmp == (int *) NULL) {
+ ffpmsg("htrans: insufficient memory");
+ return(DATA_COMPRESSION_ERR);
+ }
+ /*
+ * set up rounding and shifting masks
+ */
+ shift = 0;
+ mask = -2;
+ mask2 = mask << 1;
+ prnd = 1;
+ prnd2 = prnd << 1;
+ nrnd2 = prnd2 - 1;
+ /*
+ * do log2n reductions
+ *
+ * We're indexing a as a 2-D array with dimensions (nx,ny).
+ */
+ nxtop = nx;
+ nytop = ny;
+
+ for (k = 0; k<log2n; k++) {
+ oddx = nxtop % 2;
+ oddy = nytop % 2;
+ for (i = 0; i<nxtop-oddx; i += 2) {
+ s00 = i*ny; /* s00 is index of a[i,j] */
+ s10 = s00+ny; /* s10 is index of a[i+1,j] */
+ for (j = 0; j<nytop-oddy; j += 2) {
+ /*
+ * Divide h0,hx,hy,hc by 2 (1 the first time through).
+ */
+ h0 = (a[s10+1] + a[s10] + a[s00+1] + a[s00]) >> 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<nytop-oddy; j += 2) {
+ h0 = (a[s00+1] + a[s00]) << (1-shift);
+ hy = (a[s00+1] - a[s00]) << (1-shift);
+ a[s00+1] = ( (hy>=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<nxtop; i++) {
+ shuffle(&a[ny*i],nytop,1,tmp);
+ }
+ for (j = 0; j<nytop; j++) {
+ shuffle(&a[j],nxtop,ny,tmp);
+ }
+ /*
+ * image size reduced by 2 (round up if odd)
+ */
+ nxtop = (nxtop+1)>>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<<log2n) ) {
+ log2n += 1;
+ }
+ /*
+ * get temporary storage for shuffling elements
+ */
+ tmp = (LONGLONG *) malloc(((nmax+1)/2)*sizeof(LONGLONG));
+ if(tmp == (LONGLONG *) NULL) {
+ ffpmsg("htrans64: insufficient memory");
+ return(DATA_COMPRESSION_ERR);
+ }
+ /*
+ * set up rounding and shifting masks
+ */
+ shift = 0;
+ mask = (LONGLONG) -2;
+ mask2 = mask << 1;
+ prnd = (LONGLONG) 1;
+ prnd2 = prnd << 1;
+ nrnd2 = prnd2 - 1;
+ /*
+ * do log2n reductions
+ *
+ * We're indexing a as a 2-D array with dimensions (nx,ny).
+ */
+ nxtop = nx;
+ nytop = ny;
+
+ for (k = 0; k<log2n; k++) {
+ oddx = nxtop % 2;
+ oddy = nytop % 2;
+ for (i = 0; i<nxtop-oddx; i += 2) {
+ s00 = i*ny; /* s00 is index of a[i,j] */
+ s10 = s00+ny; /* s10 is index of a[i+1,j] */
+ for (j = 0; j<nytop-oddy; j += 2) {
+ /*
+ * Divide h0,hx,hy,hc by 2 (1 the first time through).
+ */
+ h0 = (a[s10+1] + a[s10] + a[s00+1] + a[s00]) >> 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<nytop-oddy; j += 2) {
+ h0 = (a[s00+1] + a[s00]) << (1-shift);
+ hy = (a[s00+1] - a[s00]) << (1-shift);
+ a[s00+1] = ( (hy>=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<nxtop; i++) {
+ shuffle64(&a[ny*i],nytop,1,tmp);
+ }
+ for (j = 0; j<nytop; j++) {
+ shuffle64(&a[j],nxtop,ny,tmp);
+ }
+ /*
+ * image size reduced by 2 (round up if odd)
+ */
+ nxtop = (nxtop+1)>>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; i<n; i += 2) {
+ *p1 = *p2;
+ p1 += n2;
+ p2 += (n2+n2);
+ }
+ /*
+ * put odd elements into 2nd half
+ */
+ pt = tmp;
+ for (i = 1; i<n; i += 2) {
+ *p1 = *pt;
+ p1 += n2;
+ pt += 1;
+ }
+}
+/* ######################################################################### */
+static void
+shuffle64(LONGLONG a[], int n, int n2, LONGLONG tmp[])
+{
+
+/*
+LONGLONG a[]; array to shuffle
+int n; number of elements to shuffle
+int n2; second dimension
+LONGLONG tmp[]; scratch storage
+*/
+
+int i;
+LONGLONG *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; i<n; i += 2) {
+ *p1 = *p2;
+ p1 += n2;
+ p2 += (n2+n2);
+ }
+ /*
+ * put odd elements into 2nd half
+ */
+ pt = tmp;
+ for (i = 1; i<n; i += 2) {
+ *p1 = *pt;
+ p1 += n2;
+ pt += 1;
+ }
+}
+/* ######################################################################### */
+/* ######################################################################### */
+/* Copyright (c) 1993 Association of Universities for Research
+ * in Astronomy. All rights reserved. Produced under National
+ * Aeronautics and Space Administration Contract No. NAS5-26555.
+ */
+/* digitize.c digitize H-transform
+ *
+ * Programmer: R. White Date: 11 March 1991
+ */
+
+/* ######################################################################### */
+static void
+digitize(int a[], int nx, int ny, int scale)
+{
+int d, *p;
+
+ /*
+ * round to multiple of scale
+ */
+ if (scale <= 1) return;
+ d=(scale+1)/2-1;
+ for (p=a; p <= &a[nx*ny-1]; p++) *p = ((*p>0) ? (*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;
+
+ 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) {
+ ffpmsg("encode: insufficient memory");
+ return(DATA_COMPRESSION_ERR);
+ }
+ nsign = 0;
+ bits_to_go = 8;
+ signbits[0] = 0;
+ for (i=0; i<nel; i++) {
+ if (a[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<nel; i++) {
+ q = (j>=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
+ */
+
+ /* 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] = (int) (log((float) (vmax[q]+1))/log(2.0)+0.5);
+ if ( (vmax[q]+1) > (1<<nbitplanes[q]) ) {
+ nbitplanes[q] += 1;
+ }
+ }
+*/
+
+ /*
+ * write nbitplanes
+ */
+ if (0 == qwrite(outfile, (char *) nbitplanes, sizeof(nbitplanes))) {
+ *nlength = noutchar;
+ ffpmsg("encode: output buffer too small");
+ return(DATA_COMPRESSION_ERR);
+ }
+
+ /*
+ * write coded array
+ */
+ stat = doencode(outfile, a, nx, ny, nbitplanes);
+ /*
+ * write sign bits
+ */
+
+ if (nsign > 0) {
+
+ if ( 0 == qwrite(outfile, (char *) signbits, nsign)) {
+ free(signbits);
+ *nlength = noutchar;
+ ffpmsg("encode: output buffer too small");
+ return(DATA_COMPRESSION_ERR);
+ }
+ }
+
+ free(signbits);
+ *nlength = noutchar;
+
+ if (noutchar >= noutmax) {
+ ffpmsg("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;
+
+ 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) {
+ ffpmsg("encode64: insufficient memory");
+ return(DATA_COMPRESSION_ERR);
+ }
+ nsign = 0;
+ bits_to_go = 8;
+ signbits[0] = 0;
+ for (i=0; i<nel; i++) {
+ if (a[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<nel; i++) {
+ q = (j>=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
+ */
+
+ /* 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)<<nbitplanes[q]) ) {
+ nbitplanes[q] += 1;
+ }
+ }
+*/
+
+ /*
+ * write nbitplanes
+ */
+
+ if (0 == qwrite(outfile, (char *) nbitplanes, sizeof(nbitplanes))) {
+ *nlength = noutchar;
+ ffpmsg("encode: output buffer too small");
+ return(DATA_COMPRESSION_ERR);
+ }
+
+ /*
+ * write coded array
+ */
+ stat = doencode64(outfile, a, nx, ny, nbitplanes);
+ /*
+ * write sign bits
+ */
+
+ if (nsign > 0) {
+
+ if ( 0 == qwrite(outfile, (char *) signbits, nsign)) {
+ free(signbits);
+ *nlength = noutchar;
+ ffpmsg("encode: output buffer too small");
+ return(DATA_COMPRESSION_ERR);
+ }
+ }
+
+ free(signbits);
+ *nlength = noutchar;
+
+ if (noutchar >= noutmax) {
+ ffpmsg("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 *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;
+
+ 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;
+
+ 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(void)
+{
+ 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)
+{
+ /* AND mask for the right-most n bits */
+ static int mask[9] = {0, 1, 3, 7, 15, 31, 63, 127, 255};
+ /*
+ * insert bits at end of buffer
+ */
+ buffer2 <<= n;
+/* buffer2 |= ( bits & ((1<<n)-1) ); */
+ buffer2 |= ( bits & (*(mask+n)) );
+ bits_to_go2 -= n;
+ if (bits_to_go2 <= 0) {
+ /*
+ * buffer2 full, put out top 8 bits
+ */
+
+ outfile[noutchar] = ((buffer2>>(-bits_to_go2)) & 0xff);
+
+ if (noutchar < noutmax) noutchar++;
+
+ bits_to_go2 += 8;
+ }
+ bitcount += n;
+}
+/* ######################################################################### */
+/* OUTPUT a 4 bit nybble */
+static void
+output_nybble(char *outfile, int bits)
+{
+ /*
+ * insert 4 bits at end of buffer
+ */
+ buffer2 = (buffer2<<4) | ( bits & 15 );
+ bits_to_go2 -= 4;
+ if (bits_to_go2 <= 0) {
+ /*
+ * buffer2 full, put out top 8 bits
+ */
+
+ outfile[noutchar] = ((buffer2>>(-bits_to_go2)) & 0xff);
+
+ if (noutchar < noutmax) noutchar++;
+
+ bits_to_go2 += 8;
+ }
+ bitcount += 4;
+}
+/* ############################################################################ */
+/* OUTPUT array of 4 BITS */
+
+static void output_nnybble(char *outfile, int n, unsigned char array[])
+{
+ /* pack the 4 lower bits in each element of the array into the outfile array */
+
+int ii, jj, kk = 0, shift;
+
+ if (n == 1) {
+ output_nybble(outfile, (int) array[0]);
+ return;
+ }
+/* forcing byte alignment doesn;t help, and even makes it go slightly slower
+if (bits_to_go2 != 8)
+ output_nbits(outfile, kk, bits_to_go2);
+*/
+ if (bits_to_go2 <= 4)
+ {
+ /* just room for 1 nybble; write it out separately */
+ output_nybble(outfile, array[0]);
+ kk++; /* index to next array element */
+
+ if (n == 2) /* only 1 more nybble to write out */
+ {
+ output_nybble(outfile, (int) array[1]);
+ return;
+ }
+ }
+
+
+ /* bits_to_go2 is now in the range 5 - 8 */
+ shift = 8 - bits_to_go2;
+
+ /* now write out pairs of nybbles; this does not affect value of bits_to_go2 */
+ jj = (n - kk) / 2;
+
+ if (bits_to_go2 == 8) {
+ /* special case if nybbles are aligned on byte boundary */
+ /* this actually seems to make very little differnece in speed */
+ buffer2 = 0;
+ for (ii = 0; ii < jj; ii++)
+ {
+ outfile[noutchar] = ((array[kk] & 15)<<4) | (array[kk+1] & 15);
+ kk += 2;
+ noutchar++;
+ }
+ } else {
+ for (ii = 0; ii < jj; ii++)
+ {
+ buffer2 = (buffer2<<8) | ((array[kk] & 15)<<4) | (array[kk+1] & 15);
+ kk += 2;
+
+ /*
+ buffer2 full, put out top 8 bits
+ */
+
+ outfile[noutchar] = ((buffer2>>shift) & 0xff);
+ noutchar++;
+ }
+ }
+
+ bitcount += (8 * (ii - 1));
+
+ /* write out last odd nybble, if present */
+ if (kk != n) output_nybble(outfile, (int) array[n - 1]);
+
+ return;
+}
+
+
+/* ######################################################################### */
+/* FLUSH OUT THE LAST BITS */
+
+static void
+done_outputing_bits(char *outfile)
+{
+ if(bits_to_go2 < 8) {
+/* putc(buffer2<<bits_to_go2,outfile); */
+
+ outfile[noutchar] = (buffer2<<bits_to_go2);
+ if (noutchar < noutmax) noutchar++;
+
+ /* count the garbage bits too */
+ bitcount += bits_to_go2;
+ }
+}
+/* ######################################################################### */
+/* ######################################################################### */
+/* Copyright (c) 1993 Association of Universities for Research
+ * in Astronomy. All rights reserved. Produced under National
+ * Aeronautics and Space Administration Contract No. NAS5-26555.
+ */
+/* qtree_encode.c Encode values in quadrant of 2-D array using binary
+ * quadtree coding for each bit plane. Assumes array is
+ * positive.
+ *
+ * Programmer: R. White Date: 15 May 1991
+ */
+
+/*
+ * Huffman code values and number of bits in each code
+ */
+static int code[16] =
+ {
+ 0x3e, 0x00, 0x01, 0x08, 0x02, 0x09, 0x1a, 0x1b,
+ 0x03, 0x1c, 0x0a, 0x1d, 0x0b, 0x1e, 0x3f, 0x0c
+ };
+static int ncode[16] =
+ {
+ 6, 3, 3, 4, 3, 4, 5, 5,
+ 3, 5, 4, 5, 4, 5, 6, 4
+ };
+
+/*
+ * variables for bit output to buffer when Huffman coding
+ */
+static int bitbuffer, bits_to_go3;
+
+/*
+ * macros to write out 4-bit nybble, Huffman code for this value
+ */
+
+
+/* ######################################################################### */
+static int
+qtree_encode(char *outfile, int a[], int n, int nqx, int nqy, int nbitplanes)
+{
+
+/*
+int 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, bmax, nqmax, nqx2, nqy2, nx, ny;
+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<<log2n)) {
+ log2n += 1;
+ }
+ /*
+ * initialize buffer point, max buffer size
+ */
+ nqx2 = (nqx+1)/2;
+ nqy2 = (nqy+1)/2;
+ bmax = (nqx2*nqy2+1)/2;
+ /*
+ * We're indexing A as a 2-D array with dimensions (nqx,nqy).
+ * Scratch is 2-D with dimensions (nqx/2,nqy/2) rounded up.
+ * Buffer is used to store string of codes for output.
+ */
+ scratch = (unsigned char *) malloc(2*bmax);
+ buffer = (unsigned char *) malloc(bmax);
+ if ((scratch == (unsigned char *) NULL) ||
+ (buffer == (unsigned char *) NULL)) {
+ ffpmsg("qtree_encode: insufficient memory");
+ return(DATA_COMPRESSION_ERR);
+ }
+ /*
+ * now encode each bit plane, starting with the top
+ */
+ for (bit=nbitplanes-1; bit >= 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<log2n; k++) {
+ qtree_reduce(scratch,ny,nx,ny,scratch);
+ nx = (nx+1)>>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<<bits_to_go3)-1),
+ bits_to_go3);
+ } else {
+ /*
+ * have to write a zero nybble if there are no 1's in array
+ */
+ output_huffman(outfile,0);
+ }
+ } else {
+ if (bits_to_go3>0) {
+ /*
+ * put out the last few bits
+ */
+ output_nbits(outfile, bitbuffer & ((1<<bits_to_go3)-1),
+ bits_to_go3);
+ }
+ for (i=b-1; i>=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<<log2n)) {
+ log2n += 1;
+ }
+ /*
+ * initialize buffer point, max buffer size
+ */
+ nqx2 = (nqx+1)/2;
+ nqy2 = (nqy+1)/2;
+ bmax = (( nqx2)* ( nqy2)+1)/2;
+ /*
+ * We're indexing A as a 2-D array with dimensions (nqx,nqy).
+ * Scratch is 2-D with dimensions (nqx/2,nqy/2) rounded up.
+ * Buffer is used to store string of codes for output.
+ */
+ scratch = (unsigned char *) malloc(2*bmax);
+ buffer = (unsigned char *) malloc(bmax);
+ if ((scratch == (unsigned char *) NULL) ||
+ (buffer == (unsigned char *) NULL)) {
+ ffpmsg("qtree_encode64: insufficient memory");
+ return(DATA_COMPRESSION_ERR);
+ }
+ /*
+ * now encode each bit plane, starting with the top
+ */
+ for (bit=nbitplanes-1; bit >= 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<log2n; k++) {
+ qtree_reduce(scratch,ny,nx,ny,scratch);
+ nx = (nx+1)>>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<<bits_to_go3)-1),
+ bits_to_go3);
+ } else {
+ /*
+ * have to write a zero nybble if there are no 1's in array
+ */
+ output_huffman(outfile,0);
+ }
+ } else {
+ if (bits_to_go3>0) {
+ /*
+ * put out the last few bits
+ */
+ output_nbits(outfile, bitbuffer & ((1<<bits_to_go3)-1),
+ bits_to_go3);
+ }
+ for (i=b-1; i>=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;
+ b1 = b0<<1;
+ b2 = b0<<2;
+ b3 = b0<<3;
+ k = 0; /* k is index of b[i/2,j/2] */
+ for (i = 0; i<nx-1; i += 2) {
+ s00 = n*i; /* s00 is index of a[i,j] */
+/* tried using s00+n directly in the statements, but this had no effect on performance */
+ s10 = s00+n; /* s10 is index of a[i+1,j] */
+ for (j = 0; j<ny-1; j += 2) {
+
+/*
+ this was not any faster..
+
+ b[k] = (a[s00] & b0) ?
+ (a[s00+1] & b0) ?
+ (a[s10] & b0) ?
+ (a[s10+1] & b0) ? 15 : 14
+ : (a[s10+1] & b0) ? 13 : 12
+ : (a[s10] & b0) ?
+ (a[s10+1] & b0) ? 11 : 10
+ : (a[s10+1] & b0) ? 9 : 8
+ : (a[s00+1] & b0) ?
+ (a[s10] & b0) ?
+ (a[s10+1] & b0) ? 7 : 6
+ : (a[s10+1] & b0) ? 5 : 4
+
+ : (a[s10] & b0) ?
+ (a[s10+1] & b0) ? 3 : 2
+ : (a[s10+1] & b0) ? 1 : 0;
+*/
+
+/*
+this alternative way of calculating b[k] was slowwer than the original code
+ if ( a[s00] & b0)
+ if ( a[s00+1] & b0)
+ if ( a[s10] & b0)
+ if ( a[s10+1] & b0)
+ b[k] = 15;
+ else
+ b[k] = 14;
+ else
+ if ( a[s10+1] & b0)
+ b[k] = 13;
+ else
+ b[k] = 12;
+ else
+ if ( a[s10] & b0)
+ if ( a[s10+1] & b0)
+ b[k] = 11;
+ else
+ b[k] = 10;
+ else
+ if ( a[s10+1] & b0)
+ b[k] = 9;
+ else
+ b[k] = 8;
+ else
+ if ( a[s00+1] & b0)
+ if ( a[s10] & b0)
+ if ( a[s10+1] & b0)
+ b[k] = 7;
+ else
+ b[k] = 6;
+ else
+ if ( a[s10+1] & b0)
+ b[k] = 5;
+ else
+ b[k] = 4;
+ else
+ if ( a[s10] & b0)
+ if ( a[s10+1] & b0)
+ b[k] = 3;
+ else
+ b[k] = 2;
+ else
+ if ( a[s10+1] & b0)
+ b[k] = 1;
+ else
+ b[k] = 0;
+*/
+
+
+
+ b[k] = ( ( a[s10+1] & b0)
+ | ((a[s10 ]<<1) & b1)
+ | ((a[s00+1]<<2) & b2)
+ | ((a[s00 ]<<3) & b3) ) >> 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<ny-1; j += 2) {
+ b[k] = ( ((a[s00+1]<<2) & b2)
+ | ((a[s00 ]<<3) & b3) ) >> 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;
+ b1 = b0<<1;
+ b2 = b0<<2;
+ b3 = b0<<3;
+ k = 0; /* k is index of b[i/2,j/2] */
+ for (i = 0; i<nx-1; i += 2) {
+ s00 = n*i; /* s00 is index of a[i,j] */
+ s10 = s00+n; /* s10 is index of a[i+1,j] */
+ for (j = 0; j<ny-1; j += 2) {
+ b[k] = (unsigned char) (( ( a[s10+1] & b0)
+ | ((a[s10 ]<<1) & b1)
+ | ((a[s00+1]<<2) & b2)
+ | ((a[s00 ]<<3) & b3) ) >> 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<ny-1; j += 2) {
+ b[k] = (unsigned char) (( ((a[s00+1]<<2) & b2)
+ | ((a[s00 ]<<3) & b3) ) >> 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<nx-1; i += 2) {
+ s00 = n*i; /* s00 is index of a[i,j] */
+ s10 = s00+n; /* s10 is index of a[i+1,j] */
+ for (j = 0; j<ny-1; j += 2) {
+ b[k] = (a[s10+1] != 0)
+ | ( (a[s10 ] != 0) << 1)
+ | ( (a[s00+1] != 0) << 2)
+ | ( (a[s00 ] != 0) << 3);
+ 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 ] != 0) << 1)
+ | ( (a[s00 ] != 0) << 3);
+ 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<ny-1; j += 2) {
+ b[k] = ( (a[s00+1] != 0) << 2)
+ | ( (a[s00 ] != 0) << 3);
+ 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 ] != 0) << 3);
+ k += 1;
+ }
+ }
+}
+
+/* ######################################################################### */
+static void
+write_bdirect(char *outfile, int a[], int n,int nqx, int nqy, unsigned char scratch[], int bit)
+{
+
+ /*
+ * Write the direct bitmap warning code
+ */
+ output_nybble(outfile,0x0);
+ /*
+ * Copy A to scratch array (again!), packing 4 bits/nybble
+ */
+ qtree_onebit(a,n,nqx,nqy,scratch,bit);
+ /*
+ * write to outfile
+ */
+/*
+int i;
+ for (i = 0; i < ((nqx+1)/2) * ((nqy+1)/2); i++) {
+ output_nybble(outfile,scratch[i]);
+ }
+*/
+ output_nnybble(outfile, ((nqx+1)/2) * ((nqy+1)/2), scratch);
+
+}
+/* ######################################################################### */
+static void
+write_bdirect64(char *outfile, LONGLONG a[], int n,int nqx, int nqy, unsigned char scratch[], int bit)
+{
+
+ /*
+ * Write the direct bitmap warning code
+ */
+ output_nybble(outfile,0x0);
+ /*
+ * Copy A to scratch array (again!), packing 4 bits/nybble
+ */
+ qtree_onebit64(a,n,nqx,nqy,scratch,bit);
+ /*
+ * write to outfile
+ */
+/*
+int i;
+ for (i = 0; i < ((nqx+1)/2) * ((nqy+1)/2); i++) {
+ output_nybble(outfile,scratch[i]);
+ }
+*/
+ output_nnybble(outfile, ((nqx+1)/2) * ((nqy+1)/2), scratch);
+}
generated by cgit (git 2.34.1) at 2025-09-20 19:17:34 -0400