/*  #########################################################################
These routines to apply the H-compress decompression 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 
  hinv.c 
  hsmooth.c 
  undigitize.c 
  decode.c 
  dodecode.c 
  qtree_decode.c 
  qread.c 
  bit_input.c


The following modifications have been made to the original code:

  - commented out redundant "include" statements
  - added the nextchar global variable 
  - changed all the 'extern' declarations to 'static', since all the routines are in
    the same source file
  - changed the first parameter in decode (and in lower level routines from a file stream
    to a char array
  - modified the myread routine, and lower level byte reading routines,  to copy 
    the input bytes to a char array, instead of reading them from a file stream
  - changed the function declarations to the more modern ANSI C style
  - changed calls to printf and perror to call the CFITSIO ffpmsg routine
  - replace "exit" statements with "return" statements

 ############################################################################  */
 
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
/*#include "fitsio2.h"*/
#define LONGLONG long long
#define DATA_DECOMPRESSION_ERR 0

static void ffpmsg(const char* str) {}

/* WDP added test to see if min and max are already defined */
#ifndef min
#define min(a,b)        (((a)<(b))?(a):(b))
#endif
#ifndef max
#define max(a,b)        (((a)>(b))?(a):(b))
#endif

static long nextchar;

static int decode(unsigned char *infile, int *a, int *nx, int *ny, int *scale);
static int decode64(unsigned char *infile, LONGLONG *a, int *nx, int *ny, int *scale);
static int hinv(int a[], int nx, int ny, int smooth ,int scale);
static int hinv64(LONGLONG a[], int nx, int ny, int smooth ,int scale);
static void undigitize(int a[], int nx, int ny, int scale);
static void undigitize64(LONGLONG a[], int nx, int ny, int scale);
static void unshuffle(int a[], int n, int n2, int tmp[]);
static void unshuffle64(LONGLONG a[], int n, int n2, LONGLONG tmp[]);
static void hsmooth(int a[], int nxtop, int nytop, int ny, int scale);
static void hsmooth64(LONGLONG a[], int nxtop, int nytop, int ny, int scale);
static void qread(unsigned char *infile,char *a, int n);
static int  readint(unsigned char *infile);
static LONGLONG readlonglong(unsigned char *infile);
static int dodecode(unsigned char *infile, int a[], int nx, int ny, unsigned char nbitplanes[3]);
static int dodecode64(unsigned char *infile, LONGLONG a[], int nx, int ny, unsigned char nbitplanes[3]);
static int qtree_decode(unsigned char *infile, int a[], int n, int nqx, int nqy, int nbitplanes);
static int qtree_decode64(unsigned char *infile, LONGLONG a[], int n, int nqx, int nqy, int nbitplanes);
static void start_inputing_bits();
static int input_bit(unsigned char *infile);
static int input_nbits(unsigned char *infile, int n);
/*  make input_nybble a separate routine, for added effiency */
/* #define input_nybble(infile)	input_nbits(infile,4) */
static int input_nybble(unsigned char *infile);
static int input_nnybble(unsigned char *infile, int n, unsigned char *array);

static void qtree_expand(unsigned char *infile, unsigned char a[], int nx, int ny, unsigned char b[]);
static void qtree_bitins(unsigned char a[], int nx, int ny, int b[], int n, int bit);
static void qtree_bitins64(unsigned char a[], int nx, int ny, LONGLONG b[], int n, int bit);
static void qtree_copy(unsigned char a[], int nx, int ny, unsigned char b[], int n);
static void read_bdirect(unsigned char *infile, int a[], int n, int nqx, int nqy, unsigned char scratch[], int bit);
static void read_bdirect64(unsigned char *infile, LONGLONG a[], int n, int nqx, int nqy, unsigned char scratch[], int bit);
static int  input_huffman(unsigned char *infile);

/* ---------------------------------------------------------------------- */
int fits_hdecompress(unsigned char *input, int smooth, int *a, int *ny, int *nx, 
                     int *scale, int *status)
{
  /* 
     decompress the input byte stream using the H-compress algorithm
  
   input  - input array of compressed bytes
   a - pre-allocated array to hold the output uncompressed image
   nx - returned X axis size
   ny - returned Y axis size

 NOTE: the nx and ny dimensions as defined within this code are reversed from
 the usual FITS notation.  ny is the fastest varying dimension, which is
 usually considered the X axis in the FITS image display

  */
int stat;

  if (*status > 0) return(*status);

	/* decode the input array */

	stat = decode(input, a, nx, ny, scale);
        *status = stat;
	if (stat) return(*status);
	
	/*
	 * Un-Digitize
	 */
	undigitize(a, *nx, *ny, *scale);

	/*
	 * Inverse H-transform
	 */
	stat = hinv(a, *nx, *ny, smooth, *scale);
        *status = stat;
	
  return(*status);
}
/* ---------------------------------------------------------------------- */
int fits_hdecompress64(unsigned char *input, int smooth, LONGLONG *a, int *ny, int *nx, 
                     int *scale, int *status)
{
  /* 
     decompress the input byte stream using the H-compress algorithm
  
   input  - input array of compressed bytes
   a - pre-allocated array to hold the output uncompressed image
   nx - returned X axis size
   ny - returned Y axis size

 NOTE: the nx and ny dimensions as defined within this code are reversed from
 the usual FITS notation.  ny is the fastest varying dimension, which is
 usually considered the X axis in the FITS image display

  */
  int stat, *iarray, ii, nval;

  if (*status > 0) return(*status);

	/* decode the input array */

	stat = decode64(input, a, nx, ny, scale);
        *status = stat;
	if (stat) return(*status);
	
	/*
	 * Un-Digitize
	 */
	undigitize64(a, *nx, *ny, *scale);

	/*
	 * Inverse H-transform
	 */
	stat = hinv64(a, *nx, *ny, smooth, *scale);

        *status = stat;
	
         /* pack the I*8 values back into an I*4 array */
        iarray = (int *) a;
	nval = (*nx) * (*ny);

	for (ii = 0; ii < nval; ii++)
	   iarray[ii] = (int) a[ii];	

  return(*status);
}

/*  ############################################################################  */
/*  ############################################################################  */

/* Copyright (c) 1993 Association of Universities for Research 
 * in Astronomy. All rights reserved. Produced under National   
 * Aeronautics and Space Administration Contract No. NAS5-26555.
 */
/* hinv.c   Inverse H-transform of NX x NY integer image
 *
 * Programmer: R. White		Date: 23 July 1993
 */

/*  ############################################################################  */
static int 
hinv(int a[], int nx, int ny, int smooth ,int scale)
/*
int smooth;    0 for no smoothing, else smooth during inversion 
int scale;     used if smoothing is specified 
*/
{
int nmax, log2n, i, j, k;
int nxtop,nytop,nxf,nyf,c;
int oddx,oddy;
int shift, bit0, bit1, bit2, mask0, mask1, mask2,
	prnd0, prnd1, prnd2, nrnd0, nrnd1, nrnd2, lowbit0, lowbit1;
int h0, hx, hy, hc;
int s10, s00;
int *tmp;

	/*
	 * log2n is log2 of max(nx,ny) rounded up to next power of 2
	 */
	nmax = (nx>ny) ? nx : ny;
	log2n = (int) (log((float) nmax)/log(2.0)+0.5);
	if ( nmax > (1<<log2n) ) {
		log2n += 1;
	}
	/*
	 * get temporary storage for shuffling elements
	 */  
	tmp = (int *) malloc(((nmax+1)/2)*sizeof(int));
	if (tmp == (int *) NULL) {
		ffpmsg("hinv: insufficient memory");
		return(DATA_DECOMPRESSION_ERR);
	}
	/*
	 * set up masks, rounding parameters
	 */
	shift  = 1;
	bit0   = 1 << (log2n - 1);
	bit1   = bit0 << 1;
	bit2   = bit0 << 2;
	mask0  = -bit0;
	mask1  = mask0 << 1;
	mask2  = mask0 << 2;
	prnd0  = bit0 >> 1;
	prnd1  = bit1 >> 1;
	prnd2  = bit2 >> 1;
	nrnd0  = prnd0 - 1;
	nrnd1  = prnd1 - 1;
	nrnd2  = prnd2 - 1;
	/*
	 * round h0 to multiple of bit2
	 */
	a[0] = (a[0] + ((a[0] >= 0) ? prnd2 : nrnd2)) & mask2;
	/*
	 * do log2n expansions
	 *
	 * We're indexing a as a 2-D array with dimensions (nx,ny).
	 */
	nxtop = 1;
	nytop = 1;
	nxf = nx;
	nyf = ny;
	c = 1<<log2n;
	for (k = log2n-1; k>=0; k--) {
		/*
		 * this somewhat cryptic code generates the sequence
		 * ntop[k-1] = (ntop[k]+1)/2, where ntop[log2n] = n
		 */
		c = c>>1;
		nxtop = nxtop<<1;
		nytop = nytop<<1;
		if (nxf <= c) { nxtop -= 1; } else { nxf -= c; }
		if (nyf <= c) { nytop -= 1; } else { nyf -= c; }
		/*
		 * double shift and fix nrnd0 (because prnd0=0) on last pass
		 */
		if (k == 0) {
			nrnd0 = 0;
			shift = 2;
		}
		/*
		 * unshuffle in each dimension to interleave coefficients
		 */
		for (i = 0; i<nxtop; i++) {
			unshuffle(&a[ny*i],nytop,1,tmp);
		}
		for (j = 0; j<nytop; j++) {
			unshuffle(&a[j],nxtop,ny,tmp);
		}
		/*
		 * smooth by interpolating coefficients if SMOOTH != 0
		 */
		if (smooth) hsmooth(a,nxtop,nytop,ny,scale);
		oddx = nxtop % 2;
		oddy = nytop % 2;
		for (i = 0; i<nxtop-oddx; i += 2) {
			s00 = ny*i;				/* 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) {
				h0 = a[s00  ];
				hx = a[s10  ];
				hy = a[s00+1];
				hc = a[s10+1];
				/*
				 * round hx and hy to multiple of bit1, hc to multiple of bit0
				 * h0 is already a multiple of bit2
				 */
				hx = (hx + ((hx >= 0) ? prnd1 : nrnd1)) & mask1;
				hy = (hy + ((hy >= 0) ? prnd1 : nrnd1)) & mask1;
				hc = (hc + ((hc >= 0) ? prnd0 : nrnd0)) & mask0;
				/*
				 * propagate bit0 of hc to hx,hy
				 */
				lowbit0 = hc & bit0;
				hx = (hx >= 0) ? (hx - lowbit0) : (hx + lowbit0);
				hy = (hy >= 0) ? (hy - lowbit0) : (hy + lowbit0);
				/*
				 * Propagate bits 0 and 1 of hc,hx,hy to h0.
				 * This could be simplified if we assume h0>0, but then
				 * the inversion would not be lossless for images with
				 * negative pixels.
				 */
				lowbit1 = (hc ^ hx ^ hy) & bit1;
				h0 = (h0 >= 0)
					? (h0 + lowbit0 - lowbit1)
					: (h0 + ((lowbit0 == 0) ? lowbit1 : (lowbit0-lowbit1)));
				/*
				 * Divide sums by 2 (4 last time)
				 */
				a[s10+1] = (h0 + hx + hy + hc) >> shift;
				a[s10  ] = (h0 + hx - hy - hc) >> shift;
				a[s00+1] = (h0 - hx + hy - hc) >> shift;
				a[s00  ] = (h0 - hx - hy + hc) >> shift;
				s00 += 2;
				s10 += 2;
			}
			if (oddy) {
				/*
				 * do last element in row if row length is odd
				 * s00+1, s10+1 are off edge
				 */
				h0 = a[s00  ];
				hx = a[s10  ];
				hx = ((hx >= 0) ? (hx+prnd1) : (hx+nrnd1)) & mask1;
				lowbit1 = hx & bit1;
				h0 = (h0 >= 0) ? (h0 - lowbit1) : (h0 + lowbit1);
				a[s10  ] = (h0 + hx) >> shift;
				a[s00  ] = (h0 - hx) >> shift;
			}
		}
		if (oddx) {
			/*
			 * do last row if column length is odd
			 * s10, s10+1 are off edge
			 */
			s00 = ny*i;
			for (j = 0; j<nytop-oddy; j += 2) {
				h0 = a[s00  ];
				hy = a[s00+1];
				hy = ((hy >= 0) ? (hy+prnd1) : (hy+nrnd1)) & mask1;
				lowbit1 = hy & bit1;
				h0 = (h0 >= 0) ? (h0 - lowbit1) : (h0 + lowbit1);
				a[s00+1] = (h0 + hy) >> shift;
				a[s00  ] = (h0 - hy) >> shift;
				s00 += 2;
			}
			if (oddy) {
				/*
				 * do corner element if both row and column lengths are odd
				 * s00+1, s10, s10+1 are off edge
				 */
				h0 = a[s00  ];
				a[s00  ] = h0 >> shift;
			}
		}
		/*
		 * divide all the masks and rounding values by 2
		 */
		bit2 = bit1;
		bit1 = bit0;
		bit0 = bit0 >> 1;
		mask1 = mask0;
		mask0 = mask0 >> 1;
		prnd1 = prnd0;
		prnd0 = prnd0 >> 1;
		nrnd1 = nrnd0;
		nrnd0 = prnd0 - 1;
	}
	free(tmp);
	return(0);
}
/*  ############################################################################  */
static int 
hinv64(LONGLONG a[], int nx, int ny, int smooth ,int scale)
/*
int smooth;    0 for no smoothing, else smooth during inversion 
int scale;     used if smoothing is specified 
*/
{
int nmax, log2n, i, j, k;
int nxtop,nytop,nxf,nyf,c;
int oddx,oddy;
int shift;
LONGLONG mask0, mask1, mask2, prnd0, prnd1, prnd2, bit0, bit1, bit2;
LONGLONG  nrnd0, nrnd1, nrnd2, lowbit0, lowbit1;
LONGLONG h0, hx, hy, hc;
int s10, s00;
LONGLONG *tmp;

	/*
	 * log2n is log2 of max(nx,ny) rounded up to next power of 2
	 */
	nmax = (nx>ny) ? nx : ny;
	log2n = (int) (log((float) nmax)/log(2.0)+0.5);
	if ( nmax > (1<<log2n) ) {
		log2n += 1;
	}
	/*
	 * get temporary storage for shuffling elements
	 */  
	tmp = (LONGLONG *) malloc(((nmax+1)/2)*sizeof(LONGLONG));
	if (tmp == (LONGLONG *) NULL) {
		ffpmsg("hinv64: insufficient memory");
		return(DATA_DECOMPRESSION_ERR);
	}
	/*
	 * set up masks, rounding parameters
	 */
	shift  = 1;
	bit0   = ((LONGLONG) 1) << (log2n - 1);
	bit1   = bit0 << 1;
	bit2   = bit0 << 2;
	mask0  = -bit0;
	mask1  = mask0 << 1;
	mask2  = mask0 << 2;
	prnd0  = bit0 >> 1;
	prnd1  = bit1 >> 1;
	prnd2  = bit2 >> 1;
	nrnd0  = prnd0 - 1;
	nrnd1  = prnd1 - 1;
	nrnd2  = prnd2 - 1;
	/*
	 * round h0 to multiple of bit2
	 */
	a[0] = (a[0] + ((a[0] >= 0) ? prnd2 : nrnd2)) & mask2;
	/*
	 * do log2n expansions
	 *
	 * We're indexing a as a 2-D array with dimensions (nx,ny).
	 */
	nxtop = 1;
	nytop = 1;
	nxf = nx;
	nyf = ny;
	c = 1<<log2n;
	for (k = log2n-1; k>=0; k--) {
		/*
		 * this somewhat cryptic code generates the sequence
		 * ntop[k-1] = (ntop[k]+1)/2, where ntop[log2n] = n
		 */
		c = c>>1;
		nxtop = nxtop<<1;
		nytop = nytop<<1;
		if (nxf <= c) { nxtop -= 1; } else { nxf -= c; }
		if (nyf <= c) { nytop -= 1; } else { nyf -= c; }
		/*
		 * double shift and fix nrnd0 (because prnd0=0) on last pass
		 */
		if (k == 0) {
			nrnd0 = 0;
			shift = 2;
		}
		/*
		 * unshuffle in each dimension to interleave coefficients
		 */
		for (i = 0; i<nxtop; i++) {
			unshuffle64(&a[ny*i],nytop,1,tmp);
		}
		for (j = 0; j<nytop; j++) {
			unshuffle64(&a[j],nxtop,ny,tmp);
		}
		/*
		 * smooth by interpolating coefficients if SMOOTH != 0
		 */
		if (smooth) hsmooth64(a,nxtop,nytop,ny,scale);
		oddx = nxtop % 2;
		oddy = nytop % 2;
		for (i = 0; i<nxtop-oddx; i += 2) {
			s00 = ny*i;				/* 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) {
				h0 = a[s00  ];
				hx = a[s10  ];
				hy = a[s00+1];
				hc = a[s10+1];
				/*
				 * round hx and hy to multiple of bit1, hc to multiple of bit0
				 * h0 is already a multiple of bit2
				 */
				hx = (hx + ((hx >= 0) ? prnd1 : nrnd1)) & mask1;
				hy = (hy + ((hy >= 0) ? prnd1 : nrnd1)) & mask1;
				hc = (hc + ((hc >= 0) ? prnd0 : nrnd0)) & mask0;
				/*
				 * propagate bit0 of hc to hx,hy
				 */
				lowbit0 = hc & bit0;
				hx = (hx >= 0) ? (hx - lowbit0) : (hx + lowbit0);
				hy = (hy >= 0) ? (hy - lowbit0) : (hy + lowbit0);
				/*
				 * Propagate bits 0 and 1 of hc,hx,hy to h0.
				 * This could be simplified if we assume h0>0, but then
				 * the inversion would not be lossless for images with
				 * negative pixels.
				 */
				lowbit1 = (hc ^ hx ^ hy) & bit1;
				h0 = (h0 >= 0)
					? (h0 + lowbit0 - lowbit1)
					: (h0 + ((lowbit0 == 0) ? lowbit1 : (lowbit0-lowbit1)));
				/*
				 * Divide sums by 2 (4 last time)
				 */
				a[s10+1] = (h0 + hx + hy + hc) >> shift;
				a[s10  ] = (h0 + hx - hy - hc) >> shift;
				a[s00+1] = (h0 - hx + hy - hc) >> shift;
				a[s00  ] = (h0 - hx - hy + hc) >> shift;
				s00 += 2;
				s10 += 2;
			}
			if (oddy) {
				/*
				 * do last element in row if row length is odd
				 * s00+1, s10+1 are off edge
				 */
				h0 = a[s00  ];
				hx = a[s10  ];
				hx = ((hx >= 0) ? (hx+prnd1) : (hx+nrnd1)) & mask1;
				lowbit1 = hx & bit1;
				h0 = (h0 >= 0) ? (h0 - lowbit1) : (h0 + lowbit1);
				a[s10  ] = (h0 + hx) >> shift;
				a[s00  ] = (h0 - hx) >> shift;
			}
		}
		if (oddx) {
			/*
			 * do last row if column length is odd
			 * s10, s10+1 are off edge
			 */
			s00 = ny*i;
			for (j = 0; j<nytop-oddy; j += 2) {
				h0 = a[s00  ];
				hy = a[s00+1];
				hy = ((hy >= 0) ? (hy+prnd1) : (hy+nrnd1)) & mask1;
				lowbit1 = hy & bit1;
				h0 = (h0 >= 0) ? (h0 - lowbit1) : (h0 + lowbit1);
				a[s00+1] = (h0 + hy) >> shift;
				a[s00  ] = (h0 - hy) >> shift;
				s00 += 2;
			}
			if (oddy) {
				/*
				 * do corner element if both row and column lengths are odd
				 * s00+1, s10, s10+1 are off edge
				 */
				h0 = a[s00  ];
				a[s00  ] = h0 >> shift;
			}
		}
		/*
		 * divide all the masks and rounding values by 2
		 */
		bit2 = bit1;
		bit1 = bit0;
		bit0 = bit0 >> 1;
		mask1 = mask0;
		mask0 = mask0 >> 1;
		prnd1 = prnd0;
		prnd0 = prnd0 >> 1;
		nrnd1 = nrnd0;
		nrnd0 = prnd0 - 1;
	}
	free(tmp);
	return(0);
}

/*  ############################################################################  */
static void
unshuffle(int a[], int n, int n2, int tmp[])
/*
int a[];	 array to shuffle					
int n;		 number of elements to shuffle	
int n2;		 second dimension					
int tmp[];	 scratch storage					
*/
{
int i;
int nhalf;
int *p1, *p2, *pt;
 
	/*
	 * copy 2nd half of array to tmp
	 */
	nhalf = (n+1)>>1;
	pt = tmp;
	p1 = &a[n2*nhalf];				/* pointer to a[i]			*/
	for (i=nhalf; i<n; i++) {
		*pt = *p1;
		p1 += n2;
		pt += 1;
	}
	/*
	 * distribute 1st half of array to even elements
	 */
	p2 = &a[ n2*(nhalf-1) ];		/* pointer to a[i]			*/
	p1 = &a[(n2*(nhalf-1))<<1];		/* pointer to a[2*i]		*/
	for (i=nhalf-1; i >= 0; i--) {
		*p1 = *p2;
		p2 -= n2;
		p1 -= (n2+n2);
	}
	/*
	 * now distribute 2nd half of array (in tmp) to odd elements
	 */
	pt = tmp;
	p1 = &a[n2];					/* pointer to a[i]			*/
	for (i=1; i<n; i += 2) {
		*p1 = *pt;
		p1 += (n2+n2);
		pt += 1;
	}
}
/*  ############################################################################  */
static void
unshuffle64(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;
int nhalf;
LONGLONG *p1, *p2, *pt;
 
	/*
	 * copy 2nd half of array to tmp
	 */
	nhalf = (n+1)>>1;
	pt = tmp;
	p1 = &a[n2*nhalf];				/* pointer to a[i]			*/
	for (i=nhalf; i<n; i++) {
		*pt = *p1;
		p1 += n2;
		pt += 1;
	}
	/*
	 * distribute 1st half of array to even elements
	 */
	p2 = &a[ n2*(nhalf-1) ];		/* pointer to a[i]			*/
	p1 = &a[(n2*(nhalf-1))<<1];		/* pointer to a[2*i]		*/
	for (i=nhalf-1; i >= 0; i--) {
		*p1 = *p2;
		p2 -= n2;
		p1 -= (n2+n2);
	}
	/*
	 * now distribute 2nd half of array (in tmp) to odd elements
	 */
	pt = tmp;
	p1 = &a[n2];					/* pointer to a[i]			*/
	for (i=1; i<n; i += 2) {
		*p1 = *pt;
		p1 += (n2+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.
 */
/* hsmooth.c	Smooth H-transform image by adjusting coefficients toward
 *				interpolated values
 *
 * Programmer: R. White		Date: 13 April 1992
 */

/*  ############################################################################  */
static void 
hsmooth(int a[], int nxtop, int nytop, int ny, int scale)
/*
int a[];			 array of H-transform coefficients		
int nxtop,nytop;	 size of coefficient block to use			
int ny;				 actual 1st dimension of array			
int scale;			 truncation scale factor that was used	
*/
{
int i, j;
int ny2, s10, s00, diff, dmax, dmin, s, smax;
int hm, h0, hp, hmm, hpm, hmp, hpp, hx2, hy2;
int m1,m2;

	/*
	 * Maximum change in coefficients is determined by scale factor.
	 * Since we rounded during division (see digitize.c), the biggest
	 * permitted change is scale/2.
	 */
	smax = (scale >> 1);
	if (smax <= 0) return;
	ny2 = ny << 1;
	/*
	 * We're indexing a as a 2-D array with dimensions (nxtop,ny) of which
	 * only (nxtop,nytop) are used.  The coefficients on the edge of the
	 * array are not adjusted (which is why the loops below start at 2
	 * instead of 0 and end at nxtop-2 instead of nxtop.)
	 */
	/*
	 * Adjust x difference hx
	 */
	for (i = 2; i<nxtop-2; i += 2) {
		s00 = ny*i;				/* s00 is index of a[i,j]	*/
		s10 = s00+ny;			/* s10 is index of a[i+1,j]	*/
		for (j = 0; j<nytop; j += 2) {
			/*
			 * hp is h0 (mean value) in next x zone, hm is h0 in previous x zone
			 */
			hm = a[s00-ny2];
			h0 = a[s00];
			hp = a[s00+ny2];
			/*
			 * diff = 8 * hx slope that would match h0 in neighboring zones
			 */
			diff = hp-hm;
			/*
			 * monotonicity constraints on diff
			 */
			dmax = max( min( (hp-h0), (h0-hm) ), 0 ) << 2;
			dmin = min( max( (hp-h0), (h0-hm) ), 0 ) << 2;
			/*
			 * if monotonicity would set slope = 0 then don't change hx.
			 * note dmax>=0, dmin<=0.
			 */
			if (dmin < dmax) {
				diff = max( min(diff, dmax), dmin);
				/*
				 * Compute change in slope limited to range +/- smax.
				 * Careful with rounding negative numbers when using
				 * shift for divide by 8.
				 */
				s = diff-(a[s10]<<3);
				s = (s>=0) ? (s>>3) : ((s+7)>>3) ;
				s = max( min(s, smax), -smax);
				a[s10] = a[s10]+s;
			}
			s00 += 2;
			s10 += 2;
		}
	}
	/*
	 * Adjust y difference hy
	 */
	for (i = 0; i<nxtop; i += 2) {
		s00 = ny*i+2;
		s10 = s00+ny;
		for (j = 2; j<nytop-2; j += 2) {
			hm = a[s00-2];
			h0 = a[s00];
			hp = a[s00+2];
			diff = hp-hm;
			dmax = max( min( (hp-h0), (h0-hm) ), 0 ) << 2;
			dmin = min( max( (hp-h0), (h0-hm) ), 0 ) << 2;
			if (dmin < dmax) {
				diff = max( min(diff, dmax), dmin);
				s = diff-(a[s00+1]<<3);
				s = (s>=0) ? (s>>3) : ((s+7)>>3) ;
				s = max( min(s, smax), -smax);
				a[s00+1] = a[s00+1]+s;
			}
			s00 += 2;
			s10 += 2;
		}
	}
	/*
	 * Adjust curvature difference hc
	 */
	for (i = 2; i<nxtop-2; i += 2) {
		s00 = ny*i+2;
		s10 = s00+ny;
		for (j = 2; j<nytop-2; j += 2) {
			/*
			 * ------------------    y
			 * | hmp |    | hpp |    |
			 * ------------------    |
			 * |     | h0 |     |    |
			 * ------------------    -------x
			 * | hmm |    | hpm |
			 * ------------------
			 */
			hmm = a[s00-ny2-2];
			hpm = a[s00+ny2-2];
			hmp = a[s00-ny2+2];
			hpp = a[s00+ny2+2];
			h0  = a[s00];
			/*
			 * diff = 64 * hc value that would match h0 in neighboring zones
			 */
			diff = hpp + hmm - hmp - hpm;
			/*
			 * 2 times x,y slopes in this zone
			 */
			hx2 = a[s10  ]<<1;
			hy2 = a[s00+1]<<1;
			/*
			 * monotonicity constraints on diff
			 */
			m1 = min(max(hpp-h0,0)-hx2-hy2, max(h0-hpm,0)+hx2-hy2);
			m2 = min(max(h0-hmp,0)-hx2+hy2, max(hmm-h0,0)+hx2+hy2);
			dmax = min(m1,m2) << 4;
			m1 = max(min(hpp-h0,0)-hx2-hy2, min(h0-hpm,0)+hx2-hy2);
			m2 = max(min(h0-hmp,0)-hx2+hy2, min(hmm-h0,0)+hx2+hy2);
			dmin = max(m1,m2) << 4;
			/*
			 * if monotonicity would set slope = 0 then don't change hc.
			 * note dmax>=0, dmin<=0.
			 */
			if (dmin < dmax) {
				diff = max( min(diff, dmax), dmin);
				/*
				 * Compute change in slope limited to range +/- smax.
				 * Careful with rounding negative numbers when using
				 * shift for divide by 64.
				 */
				s = diff-(a[s10+1]<<6);
				s = (s>=0) ? (s>>6) : ((s+63)>>6) ;
				s = max( min(s, smax), -smax);
				a[s10+1] = a[s10+1]+s;
			}
			s00 += 2;
			s10 += 2;
		}
	}
}
/*  ############################################################################  */
static void 
hsmooth64(LONGLONG a[], int nxtop, int nytop, int ny, int scale)
/*
LONGLONG a[];			 array of H-transform coefficients		
int nxtop,nytop;	 size of coefficient block to use			
int ny;				 actual 1st dimension of array			
int scale;			 truncation scale factor that was used	
*/
{
int i, j;
int ny2, s10, s00;
LONGLONG hm, h0, hp, hmm, hpm, hmp, hpp, hx2, hy2, diff, dmax, dmin, s, smax, m1, m2;

	/*
	 * Maximum change in coefficients is determined by scale factor.
	 * Since we rounded during division (see digitize.c), the biggest
	 * permitted change is scale/2.
	 */
	smax = (scale >> 1);
	if (smax <= 0) return;
	ny2 = ny << 1;
	/*
	 * We're indexing a as a 2-D array with dimensions (nxtop,ny) of which
	 * only (nxtop,nytop) are used.  The coefficients on the edge of the
	 * array are not adjusted (which is why the loops below start at 2
	 * instead of 0 and end at nxtop-2 instead of nxtop.)
	 */
	/*
	 * Adjust x difference hx
	 */
	for (i = 2; i<nxtop-2; i += 2) {
		s00 = ny*i;				/* s00 is index of a[i,j]	*/
		s10 = s00+ny;			/* s10 is index of a[i+1,j]	*/
		for (j = 0; j<nytop; j += 2) {
			/*
			 * hp is h0 (mean value) in next x zone, hm is h0 in previous x zone
			 */
			hm = a[s00-ny2];
			h0 = a[s00];
			hp = a[s00+ny2];
			/*
			 * diff = 8 * hx slope that would match h0 in neighboring zones
			 */
			diff = hp-hm;
			/*
			 * monotonicity constraints on diff
			 */
			dmax = max( min( (hp-h0), (h0-hm) ), 0 ) << 2;
			dmin = min( max( (hp-h0), (h0-hm) ), 0 ) << 2;
			/*
			 * if monotonicity would set slope = 0 then don't change hx.
			 * note dmax>=0, dmin<=0.
			 */
			if (dmin < dmax) {
				diff = max( min(diff, dmax), dmin);
				/*
				 * Compute change in slope limited to range +/- smax.
				 * Careful with rounding negative numbers when using
				 * shift for divide by 8.
				 */
				s = diff-(a[s10]<<3);
				s = (s>=0) ? (s>>3) : ((s+7)>>3) ;
				s = max( min(s, smax), -smax);
				a[s10] = a[s10]+s;
			}
			s00 += 2;
			s10 += 2;
		}
	}
	/*
	 * Adjust y difference hy
	 */
	for (i = 0; i<nxtop; i += 2) {
		s00 = ny*i+2;
		s10 = s00+ny;
		for (j = 2; j<nytop-2; j += 2) {
			hm = a[s00-2];
			h0 = a[s00];
			hp = a[s00+2];
			diff = hp-hm;
			dmax = max( min( (hp-h0), (h0-hm) ), 0 ) << 2;
			dmin = min( max( (hp-h0), (h0-hm) ), 0 ) << 2;
			if (dmin < dmax) {
				diff = max( min(diff, dmax), dmin);
				s = diff-(a[s00+1]<<3);
				s = (s>=0) ? (s>>3) : ((s+7)>>3) ;
				s = max( min(s, smax), -smax);
				a[s00+1] = a[s00+1]+s;
			}
			s00 += 2;
			s10 += 2;
		}
	}
	/*
	 * Adjust curvature difference hc
	 */
	for (i = 2; i<nxtop-2; i += 2) {
		s00 = ny*i+2;
		s10 = s00+ny;
		for (j = 2; j<nytop-2; j += 2) {
			/*
			 * ------------------    y
			 * | hmp |    | hpp |    |
			 * ------------------    |
			 * |     | h0 |     |    |
			 * ------------------    -------x
			 * | hmm |    | hpm |
			 * ------------------
			 */
			hmm = a[s00-ny2-2];
			hpm = a[s00+ny2-2];
			hmp = a[s00-ny2+2];
			hpp = a[s00+ny2+2];
			h0  = a[s00];
			/*
			 * diff = 64 * hc value that would match h0 in neighboring zones
			 */
			diff = hpp + hmm - hmp - hpm;
			/*
			 * 2 times x,y slopes in this zone
			 */
			hx2 = a[s10  ]<<1;
			hy2 = a[s00+1]<<1;
			/*
			 * monotonicity constraints on diff
			 */
			m1 = min(max(hpp-h0,0)-hx2-hy2, max(h0-hpm,0)+hx2-hy2);
			m2 = min(max(h0-hmp,0)-hx2+hy2, max(hmm-h0,0)+hx2+hy2);
			dmax = min(m1,m2) << 4;
			m1 = max(min(hpp-h0,0)-hx2-hy2, min(h0-hpm,0)+hx2-hy2);
			m2 = max(min(h0-hmp,0)-hx2+hy2, min(hmm-h0,0)+hx2+hy2);
			dmin = max(m1,m2) << 4;
			/*
			 * if monotonicity would set slope = 0 then don't change hc.
			 * note dmax>=0, dmin<=0.
			 */
			if (dmin < dmax) {
				diff = max( min(diff, dmax), dmin);
				/*
				 * Compute change in slope limited to range +/- smax.
				 * Careful with rounding negative numbers when using
				 * shift for divide by 64.
				 */
				s = diff-(a[s10+1]<<6);
				s = (s>=0) ? (s>>6) : ((s+63)>>6) ;
				s = max( min(s, smax), -smax);
				a[s10+1] = a[s10+1]+s;
			}
			s00 += 2;
			s10 += 2;
		}
	}
}


/*  ############################################################################  */
/*  ############################################################################  */
/* Copyright (c) 1993 Association of Universities for Research 
 * in Astronomy. All rights reserved. Produced under National   
 * Aeronautics and Space Administration Contract No. NAS5-26555.
 */
/* undigitize.c		undigitize H-transform
 *
 * Programmer: R. White		Date: 9 May 1991
 */

/*  ############################################################################  */
static void
undigitize(int a[], int nx, int ny, int scale)
{
int *p;

	/*
	 * multiply by scale
	 */
	if (scale <= 1) return;
	for (p=a; p <= &a[nx*ny-1]; p++) *p = (*p)*scale;
}
/*  ############################################################################  */
static void
undigitize64(LONGLONG a[], int nx, int ny, int scale)
{
LONGLONG *p, scale64;

	/*
	 * multiply by scale
	 */
	if (scale <= 1) return;
	scale64 = (LONGLONG) scale;   /* use a 64-bit int for efficiency in the big loop */
	
	for (p=a; p <= &a[nx*ny-1]; p++) *p = (*p)*scale64;
}

/*  ############################################################################  */
/*  ############################################################################  */
/* Copyright (c) 1993 Association of Universities for Research 
 * in Astronomy. All rights reserved. Produced under National   
 * Aeronautics and Space Administration Contract No. NAS5-26555.
 */
/* decode.c		read codes from infile and construct array
 *
 * Programmer: R. White		Date: 2 February 1994
 */


static char code_magic[2] = { (char)0xDD, (char)0x99 };

/*  ############################################################################  */
static int decode(unsigned char *infile, int *a, int *nx, int *ny, int *scale)
/*
char *infile;				 input file							
int  *a;				 address of output array [nx][ny]		
int  *nx,*ny;				 size of output array					
int  *scale;				 scale factor for digitization		
*/
{
LONGLONG sumall;
int nel, stat;
unsigned char nbitplanes[3];
char tmagic[2];

	/* initialize the byte read position to the beginning of the array */;
	nextchar = 0;
	
	/*
	 * File starts either with special 2-byte magic code or with
	 * FITS keyword "SIMPLE  ="
	 */
	qread(infile, tmagic, sizeof(tmagic));
	/*
	 * check for correct magic code value
	 */
	if (memcmp(tmagic,code_magic,sizeof(code_magic)) != 0) {
		ffpmsg("bad file format");
		return(DATA_DECOMPRESSION_ERR);
	}
	*nx =readint(infile);				/* x size of image			*/
	*ny =readint(infile);				/* y size of image			*/
	*scale=readint(infile);				/* scale factor for digitization	*/
	
	nel = (*nx) * (*ny);

	/* sum of all pixels	*/
	sumall=readlonglong(infile);
	/* # bits in quadrants	*/

	qread(infile, (char *) nbitplanes, sizeof(nbitplanes));

	stat = dodecode(infile, a, *nx, *ny, nbitplanes);
	/*
	 * put sum of all pixels back into pixel 0
	 */
	a[0] = (int) sumall;
	return(stat);
}
/*  ############################################################################  */
static int decode64(unsigned char *infile, LONGLONG *a, int *nx, int *ny, int *scale)
/*
char *infile;				 input file							
LONGLONG  *a;				 address of output array [nx][ny]		
int  *nx,*ny;				 size of output array					
int  *scale;				 scale factor for digitization		
*/
{
int nel, stat;
LONGLONG sumall;
unsigned char nbitplanes[3];
char tmagic[2];

	/* initialize the byte read position to the beginning of the array */;
	nextchar = 0;
	
	/*
	 * File starts either with special 2-byte magic code or with
	 * FITS keyword "SIMPLE  ="
	 */
	qread(infile, tmagic, sizeof(tmagic));
	/*
	 * check for correct magic code value
	 */
	if (memcmp(tmagic,code_magic,sizeof(code_magic)) != 0) {
		ffpmsg("bad file format");
		return(DATA_DECOMPRESSION_ERR);
	}
	*nx =readint(infile);				/* x size of image			*/
	*ny =readint(infile);				/* y size of image			*/
	*scale=readint(infile);				/* scale factor for digitization	*/
	
	nel = (*nx) * (*ny);

	/* sum of all pixels	*/
	sumall=readlonglong(infile);
	/* # bits in quadrants	*/

	qread(infile, (char *) nbitplanes, sizeof(nbitplanes));

	stat = dodecode64(infile, a, *nx, *ny, nbitplanes);
	/*
	 * put sum of all pixels back into pixel 0
	 */
	a[0] = sumall;

	return(stat);
}


/*  ############################################################################  */
/*  ############################################################################  */
/* Copyright (c) 1993 Association of Universities for Research 
 * in Astronomy. All rights reserved. Produced under National   
 * Aeronautics and Space Administration Contract No. NAS5-26555.
 */
/* dodecode.c	Decode stream of characters on infile and return array
 *
 * This version encodes the different quadrants separately
 *
 * Programmer: R. White		Date: 9 May 1991
 */

/*  ############################################################################  */
static int
dodecode(unsigned char *infile, int a[], int nx, int ny, unsigned char nbitplanes[3])

/* int a[];					 			
   int nx,ny;					 Array dimensions are [nx][ny]		
   unsigned char nbitplanes[3];		 Number of bit planes in quadrants
*/
{
int i, nel, nx2, ny2, stat;

	nel = nx*ny;
	nx2 = (nx+1)/2;
	ny2 = (ny+1)/2;

	/*
	 * initialize a to zero
	 */
	for (i=0; i<nel; i++) a[i] = 0;
	/*
	 * Initialize bit input
	 */
	start_inputing_bits();
	/*
	 * read bit planes for each quadrant
	 */
	stat = qtree_decode(infile, &a[0],          ny, nx2,  ny2,  nbitplanes[0]);
        if (stat) return(stat);
	
	stat = qtree_decode(infile, &a[ny2],        ny, nx2,  ny/2, nbitplanes[1]);
        if (stat) return(stat);
	
	stat = qtree_decode(infile, &a[ny*nx2],     ny, nx/2, ny2,  nbitplanes[1]);
        if (stat) return(stat);
	
	stat = qtree_decode(infile, &a[ny*nx2+ny2], ny, nx/2, ny/2, nbitplanes[2]);
        if (stat) return(stat);
	
	/*
	 * make sure there is an EOF symbol (nybble=0) at end
	 */
	if (input_nybble(infile) != 0) {
		ffpmsg("dodecode: bad bit plane values");
		return(DATA_DECOMPRESSION_ERR);
	}
	/*
	 * now get the sign bits
	 * Re-initialize bit input
	 */
	start_inputing_bits();
	for (i=0; i<nel; i++) {
		if (a[i]) {
			/* tried putting the input_bit code in-line here, instead of */
			/* calling the function, but it made no difference in the speed */
			if (input_bit(infile)) a[i] = -a[i];
		}
	}
	return(0);
}
/*  ############################################################################  */
static int
dodecode64(unsigned char *infile, LONGLONG a[], int nx, int ny, unsigned char nbitplanes[3])

/* LONGLONG a[];					 			
   int nx,ny;					 Array dimensions are [nx][ny]		
   unsigned char nbitplanes[3];		 Number of bit planes in quadrants
*/
{
int i, nel, nx2, ny2, stat;

	nel = nx*ny;
	nx2 = (nx+1)/2;
	ny2 = (ny+1)/2;

	/*
	 * initialize a to zero
	 */
	for (i=0; i<nel; i++) a[i] = 0;
	/*
	 * Initialize bit input
	 */
	start_inputing_bits();
	/*
	 * read bit planes for each quadrant
	 */
	stat = qtree_decode64(infile, &a[0],          ny, nx2,  ny2,  nbitplanes[0]);
        if (stat) return(stat);
	
	stat = qtree_decode64(infile, &a[ny2],        ny, nx2,  ny/2, nbitplanes[1]);
        if (stat) return(stat);
	
	stat = qtree_decode64(infile, &a[ny*nx2],     ny, nx/2, ny2,  nbitplanes[1]);
        if (stat) return(stat);
	
	stat = qtree_decode64(infile, &a[ny*nx2+ny2], ny, nx/2, ny/2, nbitplanes[2]);
        if (stat) return(stat);
	
	/*
	 * make sure there is an EOF symbol (nybble=0) at end
	 */
	if (input_nybble(infile) != 0) {
		ffpmsg("dodecode64: bad bit plane values");
		return(DATA_DECOMPRESSION_ERR);
	}
	/*
	 * now get the sign bits
	 * Re-initialize bit input
	 */
	start_inputing_bits();
	for (i=0; i<nel; i++) {
		if (a[i]) {
			if (input_bit(infile) != 0) a[i] = -a[i];
		}
	}
	return(0);
}

/*  ############################################################################  */
/*  ############################################################################  */
/* 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_decode.c	Read stream of codes from infile and construct bit planes
 *					in quadrant of 2-D array using binary quadtree coding
 *
 * Programmer: R. White		Date: 7 May 1991
 */

/*  ############################################################################  */
static int
qtree_decode(unsigned char *infile, int a[], int n, int nqx, int nqy, int nbitplanes)

/*
char *infile;
int a[];				 a is 2-D array with dimensions (n,n)	
int n;					 length of full row in a				
int nqx;				 partial length of row to decode		
int nqy;				 partial length of column (<=n)		
int nbitplanes;				 number of bitplanes to decode		
*/
{
int log2n, k, bit, b, nqmax;
int nx,ny,nfx,nfy,c;
int nqx2, nqy2;
unsigned char *scratch;

	/*
	 * log2n is log2 of max(nqx,nqy) rounded up to next power of 2
	 */
	nqmax = (nqx>nqy) ? nqx : nqy;
	log2n = (int) (log((float) nqmax)/log(2.0)+0.5);
	if (nqmax > (1<<log2n)) {
		log2n += 1;
	}
	/*
	 * allocate scratch array for working space
	 */
	nqx2=(nqx+1)/2;
	nqy2=(nqy+1)/2;
	scratch = (unsigned char *) malloc(nqx2*nqy2);
	if (scratch == (unsigned char *) NULL) {
		ffpmsg("qtree_decode: insufficient memory");
		return(DATA_DECOMPRESSION_ERR);
	}
	/*
	 * now decode each bit plane, starting at the top
	 * A is assumed to be initialized to zero
	 */
	for (bit = nbitplanes-1; bit >= 0; bit--) {
		/*
		 * Was bitplane was quadtree-coded or written directly?
		 */
		b = input_nybble(infile);

		if(b == 0) {
			/*
			 * bit map was written directly
			 */
			read_bdirect(infile,a,n,nqx,nqy,scratch,bit);
		} else if (b != 0xf) {
			ffpmsg("qtree_decode: bad format code");
			return(DATA_DECOMPRESSION_ERR);
		} else {
			/*
			 * bitmap was quadtree-coded, do log2n expansions
			 *
			 * read first code
			 */
			scratch[0] = input_huffman(infile);
			/*
			 * now do log2n expansions, reading codes from file as necessary
			 */
			nx = 1;
			ny = 1;
			nfx = nqx;
			nfy = nqy;
			c = 1<<log2n;
			for (k = 1; k<log2n; k++) {
				/*
				 * this somewhat cryptic code generates the sequence
				 * n[k-1] = (n[k]+1)/2 where n[log2n]=nqx or nqy
				 */
				c = c>>1;
				nx = nx<<1;
				ny = ny<<1;
				if (nfx <= c) { nx -= 1; } else { nfx -= c; }
				if (nfy <= c) { ny -= 1; } else { nfy -= c; }
				qtree_expand(infile,scratch,nx,ny,scratch);
			}
			/*
			 * now copy last set of 4-bit codes to bitplane bit of array a
			 */
			qtree_bitins(scratch,nqx,nqy,a,n,bit);
		}
	}
	free(scratch);
	return(0);
}
/*  ############################################################################  */
static int
qtree_decode64(unsigned char *infile, LONGLONG a[], int n, int nqx, int nqy, int nbitplanes)

/*
char *infile;
LONGLONG a[];				 a is 2-D array with dimensions (n,n)	
int n;					 length of full row in a				
int nqx;				 partial length of row to decode		
int nqy;				 partial length of column (<=n)		
int nbitplanes;				 number of bitplanes to decode		
*/
{
int log2n, k, bit, b, nqmax;
int nx,ny,nfx,nfy,c;
int nqx2, nqy2;
unsigned char *scratch;

	/*
	 * log2n is log2 of max(nqx,nqy) rounded up to next power of 2
	 */
	nqmax = (nqx>nqy) ? nqx : nqy;
	log2n = (int) (log((float) nqmax)/log(2.0)+0.5);
	if (nqmax > (1<<log2n)) {
		log2n += 1;
	}
	/*
	 * allocate scratch array for working space
	 */
	nqx2=(nqx+1)/2;
	nqy2=(nqy+1)/2;
	scratch = (unsigned char *) malloc(nqx2*nqy2);
	if (scratch == (unsigned char *) NULL) {
		ffpmsg("qtree_decode64: insufficient memory");
		return(DATA_DECOMPRESSION_ERR);
	}
	/*
	 * now decode each bit plane, starting at the top
	 * A is assumed to be initialized to zero
	 */
	for (bit = nbitplanes-1; bit >= 0; bit--) {
		/*
		 * Was bitplane was quadtree-coded or written directly?
		 */
		b = input_nybble(infile);

		if(b == 0) {
			/*
			 * bit map was written directly
			 */
			read_bdirect64(infile,a,n,nqx,nqy,scratch,bit);
		} else if (b != 0xf) {
			ffpmsg("qtree_decode64: bad format code");
			return(DATA_DECOMPRESSION_ERR);
		} else {
			/*
			 * bitmap was quadtree-coded, do log2n expansions
			 *
			 * read first code
			 */
			scratch[0] = input_huffman(infile);
			/*
			 * now do log2n expansions, reading codes from file as necessary
			 */
			nx = 1;
			ny = 1;
			nfx = nqx;
			nfy = nqy;
			c = 1<<log2n;
			for (k = 1; k<log2n; k++) {
				/*
				 * this somewhat cryptic code generates the sequence
				 * n[k-1] = (n[k]+1)/2 where n[log2n]=nqx or nqy
				 */
				c = c>>1;
				nx = nx<<1;
				ny = ny<<1;
				if (nfx <= c) { nx -= 1; } else { nfx -= c; }
				if (nfy <= c) { ny -= 1; } else { nfy -= c; }
				qtree_expand(infile,scratch,nx,ny,scratch);
			}
			/*
			 * now copy last set of 4-bit codes to bitplane bit of array a
			 */
			qtree_bitins64(scratch,nqx,nqy,a,n,bit);
		}
	}
	free(scratch);
	return(0);
}


/*  ############################################################################  */
/*
 * do one quadtree expansion step on array a[(nqx+1)/2,(nqy+1)/2]
 * results put into b[nqx,nqy] (which may be the same as a)
 */
static void
qtree_expand(unsigned char *infile, unsigned char a[], int nx, int ny, unsigned char b[])
{
int i;

	/*
	 * first copy a to b, expanding each 4-bit value
	 */
	qtree_copy(a,nx,ny,b,ny);
	/*
	 * now read new 4-bit values into b for each non-zero element
	 */
	for (i = nx*ny-1; i >= 0; i--) {
		if (b[i]) b[i] = input_huffman(infile);
	}
}

/*  ############################################################################  */
/*
 * copy 4-bit values from a[(nx+1)/2,(ny+1)/2] to b[nx,ny], expanding
 * each value to 2x2 pixels
 * a,b may be same array
 */
static void
qtree_copy(unsigned char a[], int nx, int ny, unsigned char b[], int n)
/*   int n;		declared y dimension of b */
{
int i, j, k, nx2, ny2;
int s00, s10;

	/*
	 * first copy 4-bit values to b
	 * start at end in case a,b are same array
	 */
	nx2 = (nx+1)/2;
	ny2 = (ny+1)/2;
	k = ny2*(nx2-1)+ny2-1;			/* k   is index of a[i,j]		*/
	for (i = nx2-1; i >= 0; i--) {
		s00 = 2*(n*i+ny2-1);		/* s00 is index of b[2*i,2*j]		*/
		for (j = ny2-1; j >= 0; j--) {
			b[s00] = a[k];
			k -= 1;
			s00 -= 2;
		}
	}
	/*
	 * now expand each 2x2 block
	 */
	for (i = 0; i<nx-1; i += 2) {

  /* Note:
     Unlike the case in qtree_bitins, this code runs faster on a 32-bit linux
     machine using the s10 intermediate variable, rather that using s00+n. 
     Go figure!
  */
		s00 = n*i;				/* s00 is index of b[i,j]	*/
		s10 = s00+n;				/* s10 is index of b[i+1,j]	*/

		for (j = 0; j<ny-1; j += 2) {

		    switch (b[s00]) {
		    case(0):
			b[s10+1] = 0;
			b[s10  ] = 0;
			b[s00+1] = 0;
			b[s00  ] = 0;

			break;
		    case(1):
			b[s10+1] = 1;
			b[s10  ] = 0;
			b[s00+1] = 0;
			b[s00  ] = 0;

			break;
		    case(2):
			b[s10+1] = 0;
			b[s10  ] = 1;
			b[s00+1] = 0;
			b[s00  ] = 0;

			break;
		    case(3):
			b[s10+1] = 1;
			b[s10  ] = 1;
			b[s00+1] = 0;
			b[s00  ] = 0;

			break;
		    case(4):
			b[s10+1] = 0;
			b[s10  ] = 0;
			b[s00+1] = 1;
			b[s00  ] = 0;

			break;
		    case(5):
			b[s10+1] = 1;
			b[s10  ] = 0;
			b[s00+1] = 1;
			b[s00  ] = 0;

			break;
		    case(6):
			b[s10+1] = 0;
			b[s10  ] = 1;
			b[s00+1] = 1;
			b[s00  ] = 0;

			break;
		    case(7):
			b[s10+1] = 1;
			b[s10  ] = 1;
			b[s00+1] = 1;
			b[s00  ] = 0;

			break;
		    case(8):
			b[s10+1] = 0;
			b[s10  ] = 0;
			b[s00+1] = 0;
			b[s00  ] = 1;

			break;
		    case(9):
			b[s10+1] = 1;
			b[s10  ] = 0;
			b[s00+1] = 0;
			b[s00  ] = 1;
			break;
		    case(10):
			b[s10+1] = 0;
			b[s10  ] = 1;
			b[s00+1] = 0;
			b[s00  ] = 1;

			break;
		    case(11):
			b[s10+1] = 1;
			b[s10  ] = 1;
			b[s00+1] = 0;
			b[s00  ] = 1;

			break;
		    case(12):
			b[s10+1] = 0;
			b[s10  ] = 0;
			b[s00+1] = 1;
			b[s00  ] = 1;

			break;
		    case(13):
			b[s10+1] = 1;
			b[s10  ] = 0;
			b[s00+1] = 1;
			b[s00  ] = 1;

			break;
		    case(14):
			b[s10+1] = 0;
			b[s10  ] = 1;
			b[s00+1] = 1;
			b[s00  ] = 1;

			break;
		    case(15):
			b[s10+1] = 1;
			b[s10  ] = 1;
			b[s00+1] = 1;
			b[s00  ] = 1;

			break;
		    }
/*
			b[s10+1] =  b[s00]     & 1;
			b[s10  ] = (b[s00]>>1) & 1;
			b[s00+1] = (b[s00]>>2) & 1;
			b[s00  ] = (b[s00]>>3) & 1;
*/

			s00 += 2;
			s10 += 2;
		}

		if (j < ny) {
			/*
			 * row size is odd, do last element in row
			 * s00+1, s10+1 are off edge
			 */
                        /* not worth converting this to use 16 case statements */
			b[s10  ] = (b[s00]>>1) & 1;
			b[s00  ] = (b[s00]>>3) & 1;
		}
	}
	if (i < nx) {
		/*
		 * column size is odd, do last row
		 * s10, s10+1 are off edge
		 */
		s00 = n*i;
		for (j = 0; j<ny-1; j += 2) {
                        /* not worth converting this to use 16 case statements */
			b[s00+1] = (b[s00]>>2) & 1;
			b[s00  ] = (b[s00]>>3) & 1;
			s00 += 2;
		}
		if (j < ny) {
			/*
			 * both row and column size are odd, do corner element
			 * s00+1, s10, s10+1 are off edge
			 */
                        /* not worth converting this to use 16 case statements */
			b[s00  ] = (b[s00]>>3) & 1;
		}
	}
}

/*  ############################################################################  */
/*
 * Copy 4-bit values from a[(nx+1)/2,(ny+1)/2] to b[nx,ny], expanding
 * each value to 2x2 pixels and inserting into bitplane BIT of B.
 * A,B may NOT be same array (it wouldn't make sense to be inserting
 * bits into the same array anyway.)
 */
static void
qtree_bitins(unsigned char a[], int nx, int ny, int b[], int n, int bit)
/*
   int n;		declared y dimension of b
*/
{
int i, j, k;
int s00;
int plane_val;

	plane_val = 1 << bit;
	
	/*
	 * expand each 2x2 block
	 */
	k = 0;						/* k   is index of a[i/2,j/2]	*/
	for (i = 0; i<nx-1; i += 2) {
		s00 = n*i;				/* s00 is index of b[i,j]	*/

  /* Note:
     this code appears to run very slightly faster on a 32-bit linux
     machine using s00+n rather than the s10 intermediate variable
  */
  /*		s10 = s00+n;	*/			/* s10 is index of b[i+1,j]	*/
		for (j = 0; j<ny-1; j += 2) {

		    switch (a[k]) {
		    case(0):
			break;
		    case(1):
			b[s00+n+1] |= plane_val;
			break;
		    case(2):
			b[s00+n  ] |= plane_val;
			break;
		    case(3):
			b[s00+n+1] |= plane_val;
			b[s00+n  ] |= plane_val;
			break;
		    case(4):
			b[s00+1] |= plane_val;
			break;
		    case(5):
			b[s00+n+1] |= plane_val;
			b[s00+1] |= plane_val;
			break;
		    case(6):
			b[s00+n  ] |= plane_val;
			b[s00+1] |= plane_val;
			break;
		    case(7):
			b[s00+n+1] |= plane_val;
			b[s00+n  ] |= plane_val;
			b[s00+1] |= plane_val;
			break;
		    case(8):
			b[s00  ] |= plane_val;
			break;
		    case(9):
			b[s00+n+1] |= plane_val;
			b[s00  ] |= plane_val;
			break;
		    case(10):
			b[s00+n  ] |= plane_val;
			b[s00  ] |= plane_val;
			break;
		    case(11):
			b[s00+n+1] |= plane_val;
			b[s00+n  ] |= plane_val;
			b[s00  ] |= plane_val;
			break;
		    case(12):
			b[s00+1] |= plane_val;
			b[s00  ] |= plane_val;
			break;
		    case(13):
			b[s00+n+1] |= plane_val;
			b[s00+1] |= plane_val;
			b[s00  ] |= plane_val;
			break;
		    case(14):
			b[s00+n  ] |= plane_val;
			b[s00+1] |= plane_val;
			b[s00  ] |= plane_val;
			break;
		    case(15):
			b[s00+n+1] |= plane_val;
			b[s00+n  ] |= plane_val;
			b[s00+1] |= plane_val;
			b[s00  ] |= plane_val;
			break;
		    }

/*
			b[s10+1] |= ( a[k]     & 1) << bit;
			b[s10  ] |= ((a[k]>>1) & 1) << bit;
			b[s00+1] |= ((a[k]>>2) & 1) << bit;
			b[s00  ] |= ((a[k]>>3) & 1) << bit;
*/
			s00 += 2;
/*			s10 += 2; */
			k += 1;
		}
		if (j < ny) {
			/*
			 * row size is odd, do last element in row
			 * s00+1, s10+1 are off edge
			 */

		    switch (a[k]) {
		    case(0):
			break;
		    case(1):
			break;
		    case(2):
			b[s00+n  ] |= plane_val;
			break;
		    case(3):
			b[s00+n  ] |= plane_val;
			break;
		    case(4):
			break;
		    case(5):
			break;
		    case(6):
			b[s00+n  ] |= plane_val;
			break;
		    case(7):
			b[s00+n  ] |= plane_val;
			break;
		    case(8):
			b[s00  ] |= plane_val;
			break;
		    case(9):
			b[s00  ] |= plane_val;
			break;
		    case(10):
			b[s00+n  ] |= plane_val;
			b[s00  ] |= plane_val;
			break;
		    case(11):
			b[s00+n  ] |= plane_val;
			b[s00  ] |= plane_val;
			break;
		    case(12):
			b[s00  ] |= plane_val;
			break;
		    case(13):
			b[s00  ] |= plane_val;
			break;
		    case(14):
			b[s00+n  ] |= plane_val;
			b[s00  ] |= plane_val;
			break;
		    case(15):
			b[s00+n  ] |= plane_val;
			b[s00  ] |= plane_val;
			break;
		    }

/*
			b[s10  ] |= ((a[k]>>1) & 1) << bit;
			b[s00  ] |= ((a[k]>>3) & 1) << bit;
*/
			k += 1;
		}
	}
	if (i < nx) {
		/*
		 * column size is odd, do last row
		 * s10, s10+1 are off edge
		 */
		s00 = n*i;
		for (j = 0; j<ny-1; j += 2) {

		    switch (a[k]) {
		    case(0):
			break;
		    case(1):
			break;
		    case(2):
			break;
		    case(3):
			break;
		    case(4):
			b[s00+1] |= plane_val;
			break;
		    case(5):
			b[s00+1] |= plane_val;
			break;
		    case(6):
			b[s00+1] |= plane_val;
			break;
		    case(7):
			b[s00+1] |= plane_val;
			break;
		    case(8):
			b[s00  ] |= plane_val;
			break;
		    case(9):
			b[s00  ] |= plane_val;
			break;
		    case(10):
			b[s00  ] |= plane_val;
			break;
		    case(11):
			b[s00  ] |= plane_val;
			break;
		    case(12):
			b[s00+1] |= plane_val;
			b[s00  ] |= plane_val;
			break;
		    case(13):
			b[s00+1] |= plane_val;
			b[s00  ] |= plane_val;
			break;
		    case(14):
			b[s00+1] |= plane_val;
			b[s00  ] |= plane_val;
			break;
		    case(15):
			b[s00+1] |= plane_val;
			b[s00  ] |= plane_val;
			break;
		    }

/*
			b[s00+1] |= ((a[k]>>2) & 1) << bit;
			b[s00  ] |= ((a[k]>>3) & 1) << bit;
*/

			s00 += 2;
			k += 1;
		}
		if (j < ny) {
			/*
			 * both row and column size are odd, do corner element
			 * s00+1, s10, s10+1 are off edge
			 */

		    switch (a[k]) {
		    case(0):
			break;
		    case(1):
			break;
		    case(2):
			break;
		    case(3):
			break;
		    case(4):
			break;
		    case(5):
			break;
		    case(6):
			break;
		    case(7):
			break;
		    case(8):
			b[s00  ] |= plane_val;
			break;
		    case(9):
			b[s00  ] |= plane_val;
			break;
		    case(10):
			b[s00  ] |= plane_val;
			break;
		    case(11):
			b[s00  ] |= plane_val;
			break;
		    case(12):
			b[s00  ] |= plane_val;
			break;
		    case(13):
			b[s00  ] |= plane_val;
			break;
		    case(14):
			b[s00  ] |= plane_val;
			break;
		    case(15):
			b[s00  ] |= plane_val;
			break;
		    }

/*
			b[s00  ] |= ((a[k]>>3) & 1) << bit;
*/
			k += 1;
		}
	}
}
/*  ############################################################################  */
/*
 * Copy 4-bit values from a[(nx+1)/2,(ny+1)/2] to b[nx,ny], expanding
 * each value to 2x2 pixels and inserting into bitplane BIT of B.
 * A,B may NOT be same array (it wouldn't make sense to be inserting
 * bits into the same array anyway.)
 */
static void
qtree_bitins64(unsigned char a[], int nx, int ny, LONGLONG b[], int n, int bit)
/*
   int n;		declared y dimension of b
*/
{
int i, j, k;
int s00;
int plane_val;

	plane_val = 1 << bit;

	/*
	 * expand each 2x2 block
	 */
	k = 0;							/* k   is index of a[i/2,j/2]	*/
	for (i = 0; i<nx-1; i += 2) {
		s00 = n*i;					/* s00 is index of b[i,j]		*/

  /* Note:
     this code appears to run very slightly faster on a 32-bit linux
     machine using s00+n rather than the s10 intermediate variable
  */
  /*		s10 = s00+n;	*/			/* s10 is index of b[i+1,j]	*/
		for (j = 0; j<ny-1; j += 2) {

		    switch (a[k]) {
		    case(0):
			break;
		    case(1):
			b[s00+n+1] |= plane_val;
			break;
		    case(2):
			b[s00+n  ] |= plane_val;
			break;
		    case(3):
			b[s00+n+1] |= plane_val;
			b[s00+n  ] |= plane_val;
			break;
		    case(4):
			b[s00+1] |= plane_val;
			break;
		    case(5):
			b[s00+n+1] |= plane_val;
			b[s00+1] |= plane_val;
			break;
		    case(6):
			b[s00+n  ] |= plane_val;
			b[s00+1] |= plane_val;
			break;
		    case(7):
			b[s00+n+1] |= plane_val;
			b[s00+n  ] |= plane_val;
			b[s00+1] |= plane_val;
			break;
		    case(8):
			b[s00  ] |= plane_val;
			break;
		    case(9):
			b[s00+n+1] |= plane_val;
			b[s00  ] |= plane_val;
			break;
		    case(10):
			b[s00+n  ] |= plane_val;
			b[s00  ] |= plane_val;
			break;
		    case(11):
			b[s00+n+1] |= plane_val;
			b[s00+n  ] |= plane_val;
			b[s00  ] |= plane_val;
			break;
		    case(12):
			b[s00+1] |= plane_val;
			b[s00  ] |= plane_val;
			break;
		    case(13):
			b[s00+n+1] |= plane_val;
			b[s00+1] |= plane_val;
			b[s00  ] |= plane_val;
			break;
		    case(14):
			b[s00+n  ] |= plane_val;
			b[s00+1] |= plane_val;
			b[s00  ] |= plane_val;
			break;
		    case(15):
			b[s00+n+1] |= plane_val;
			b[s00+n  ] |= plane_val;
			b[s00+1] |= plane_val;
			b[s00  ] |= plane_val;
			break;
		    }

/*
			b[s10+1] |= ((LONGLONG) ( a[k]     & 1)) << bit;
			b[s10  ] |= ((((LONGLONG)a[k])>>1) & 1) << bit;
			b[s00+1] |= ((((LONGLONG)a[k])>>2) & 1) << bit;
			b[s00  ] |= ((((LONGLONG)a[k])>>3) & 1) << bit;
*/
			s00 += 2;
/*			s10 += 2;  */
			k += 1;
		}
		if (j < ny) {
			/*
			 * row size is odd, do last element in row
			 * s00+1, s10+1 are off edge
			 */

		    switch (a[k]) {
		    case(0):
			break;
		    case(1):
			break;
		    case(2):
			b[s00+n  ] |= plane_val;
			break;
		    case(3):
			b[s00+n  ] |= plane_val;
			break;
		    case(4):
			break;
		    case(5):
			break;
		    case(6):
			b[s00+n  ] |= plane_val;
			break;
		    case(7):
			b[s00+n  ] |= plane_val;
			break;
		    case(8):
			b[s00  ] |= plane_val;
			break;
		    case(9):
			b[s00  ] |= plane_val;
			break;
		    case(10):
			b[s00+n  ] |= plane_val;
			b[s00  ] |= plane_val;
			break;
		    case(11):
			b[s00+n  ] |= plane_val;
			b[s00  ] |= plane_val;
			break;
		    case(12):
			b[s00  ] |= plane_val;
			break;
		    case(13):
			b[s00  ] |= plane_val;
			break;
		    case(14):
			b[s00+n  ] |= plane_val;
			b[s00  ] |= plane_val;
			break;
		    case(15):
			b[s00+n  ] |= plane_val;
			b[s00  ] |= plane_val;
			break;
		    }
/*
			b[s10  ] |= ((((LONGLONG)a[k])>>1) & 1) << bit;
			b[s00  ] |= ((((LONGLONG)a[k])>>3) & 1) << bit;
*/
			k += 1;
		}
	}
	if (i < nx) {
		/*
		 * column size is odd, do last row
		 * s10, s10+1 are off edge
		 */
		s00 = n*i;
		for (j = 0; j<ny-1; j += 2) {

		    switch (a[k]) {
		    case(0):
			break;
		    case(1):
			break;
		    case(2):
			break;
		    case(3):
			break;
		    case(4):
			b[s00+1] |= plane_val;
			break;
		    case(5):
			b[s00+1] |= plane_val;
			break;
		    case(6):
			b[s00+1] |= plane_val;
			break;
		    case(7):
			b[s00+1] |= plane_val;
			break;
		    case(8):
			b[s00  ] |= plane_val;
			break;
		    case(9):
			b[s00  ] |= plane_val;
			break;
		    case(10):
			b[s00  ] |= plane_val;
			break;
		    case(11):
			b[s00  ] |= plane_val;
			break;
		    case(12):
			b[s00+1] |= plane_val;
			b[s00  ] |= plane_val;
			break;
		    case(13):
			b[s00+1] |= plane_val;
			b[s00  ] |= plane_val;
			break;
		    case(14):
			b[s00+1] |= plane_val;
			b[s00  ] |= plane_val;
			break;
		    case(15):
			b[s00+1] |= plane_val;
			b[s00  ] |= plane_val;
			break;
		    }

/*
			b[s00+1] |= ((((LONGLONG)a[k])>>2) & 1) << bit;
			b[s00  ] |= ((((LONGLONG)a[k])>>3) & 1) << bit;
*/
			s00 += 2;
			k += 1;
		}
		if (j < ny) {
			/*
			 * both row and column size are odd, do corner element
			 * s00+1, s10, s10+1 are off edge
			 */

		    switch (a[k]) {
		    case(0):
			break;
		    case(1):
			break;
		    case(2):
			break;
		    case(3):
			break;
		    case(4):
			break;
		    case(5):
			break;
		    case(6):
			break;
		    case(7):
			break;
		    case(8):
			b[s00  ] |= plane_val;
			break;
		    case(9):
			b[s00  ] |= plane_val;
			break;
		    case(10):
			b[s00  ] |= plane_val;
			break;
		    case(11):
			b[s00  ] |= plane_val;
			break;
		    case(12):
			b[s00  ] |= plane_val;
			break;
		    case(13):
			b[s00  ] |= plane_val;
			break;
		    case(14):
			b[s00  ] |= plane_val;
			break;
		    case(15):
			b[s00  ] |= plane_val;
			break;
		    }
/*
			b[s00  ] |= ((((LONGLONG)a[k])>>3) & 1) << bit;
*/
			k += 1;
		}
	}
}

/*  ############################################################################  */
static void
read_bdirect(unsigned char *infile, int a[], int n, int nqx, int nqy, unsigned char scratch[], int bit)
{


	/*
	 * read bit image packed 4 pixels/nybble
	 */
/*
        int i;
	for (i = 0; i < ((nqx+1)/2) * ((nqy+1)/2); i++) {
		scratch[i] = input_nybble(infile);
	}
*/
        input_nnybble(infile, ((nqx+1)/2) * ((nqy+1)/2), scratch);
	
	/*
	 * insert in bitplane BIT of image A
	 */
	qtree_bitins(scratch,nqx,nqy,a,n,bit);
}
/*  ############################################################################  */
static void
read_bdirect64(unsigned char *infile, LONGLONG a[], int n, int nqx, int nqy, unsigned char scratch[], int bit)
{

	/*
	 * read bit image packed 4 pixels/nybble
	 */
/*
        int i;
	for (i = 0; i < ((nqx+1)/2) * ((nqy+1)/2); i++) {
		scratch[i] = input_nybble(infile);
	}
*/
        input_nnybble(infile, ((nqx+1)/2) * ((nqy+1)/2), scratch);

	/*
	 * insert in bitplane BIT of image A
	 */
	qtree_bitins64(scratch,nqx,nqy,a,n,bit);
}

/*  ############################################################################  */
/*
 * Huffman decoding for fixed codes
 *
 * Coded values range from 0-15
 *
 * Huffman code values (hex):
 *
 *	3e, 00, 01, 08, 02, 09, 1a, 1b,
 *	03, 1c, 0a, 1d, 0b, 1e, 3f, 0c
 *
 * and number of bits in each code:
 *
 *	6,  3,  3,  4,  3,  4,  5,  5,
 *	3,  5,  4,  5,  4,  5,  6,  4
 */
static int input_huffman(unsigned char *infile)
{
int c;

	/*
	 * get first 3 bits to start
	 */
	c = input_nbits(infile,3);
	if (c < 4) {
		/*
		 * this is all we need
		 * return 1,2,4,8 for c=0,1,2,3
		 */
		return(1<<c);
	}
	/*
	 * get the next bit
	 */
	c = input_bit(infile) | (c<<1);
	if (c < 13) {
		/*
		 * OK, 4 bits is enough
		 */
		switch (c) {
			case  8 : return(3);
			case  9 : return(5);
			case 10 : return(10);
			case 11 : return(12);
			case 12 : return(15);
		}
	}
	/*
	 * get yet another bit
	 */
	c = input_bit(infile) | (c<<1);
	if (c < 31) {
		/*
		 * OK, 5 bits is enough
		 */
		switch (c) {
			case 26 : return(6);
			case 27 : return(7);
			case 28 : return(9);
			case 29 : return(11);
			case 30 : return(13);
		}
	}
	/*
	 * need the 6th bit
	 */
	c = input_bit(infile) | (c<<1);
	if (c == 62) {
		return(0);
	} else {
		return(14);
	}
}

/*  ############################################################################  */
/*  ############################################################################  */
/* Copyright (c) 1993 Association of Universities for Research 
 * in Astronomy. All rights reserved. Produced under National   
 * Aeronautics and Space Administration Contract No. NAS5-26555.
 */
/* qread.c	Read binary data
 *
 * Programmer: R. White		Date: 11 March 1991
 */

static int readint(unsigned char *infile)
{
int a,i;
unsigned char b[4];

	/* Read integer A one byte at a time from infile.
	 *
	 * This is portable from Vax to Sun since it eliminates the
	 * need for byte-swapping.
	 *
         *  This routine is only called to read the first 3 values
	 *  in the compressed file, so it doesn't have to be 
	 *  super-efficient
	 */
	for (i=0; i<4; i++) qread(infile,(char *) &b[i],1);
	a = b[0];
	for (i=1; i<4; i++) a = (a<<8) + b[i];
	return(a);
}

/*  ############################################################################  */
static LONGLONG readlonglong(unsigned char *infile)
{
int i;
LONGLONG a;
unsigned char b[8];

	/* Read integer A one byte at a time from infile.
	 *
	 * This is portable from Vax to Sun since it eliminates the
	 * need for byte-swapping.
	 *
         *  This routine is only called to read the first 3 values
	 *  in the compressed file, so it doesn't have to be 
	 *  super-efficient
	 */
	for (i=0; i<8; i++) qread(infile,(char *) &b[i],1);
	a = b[0];
	for (i=1; i<8; i++) a = (a<<8) + b[i];
	return(a);
}

/*  ############################################################################  */
static void qread(unsigned char *file, char buffer[], int n)
{
    /*
     * read n bytes from file into buffer
     *
     */

    memcpy(buffer, &file[nextchar], n);
    nextchar += n;
}

/*  ############################################################################  */
/*  ############################################################################  */
/* Copyright (c) 1993 Association of Universities for Research
 * in Astronomy. All rights reserved. Produced under National
 * Aeronautics and Space Administration Contract No. NAS5-26555.
 */

/* BIT INPUT ROUTINES */

/* THE BIT BUFFER */

static int buffer2;			/* Bits waiting to be input	*/
static int bits_to_go;			/* Number of bits still in buffer */

/* INITIALIZE BIT INPUT */

/*  ############################################################################  */
static void start_inputing_bits()
{
	/*
	 * Buffer starts out with no bits in it
	 */
	bits_to_go = 0;
}

/*  ############################################################################  */
/* INPUT A BIT */

static int input_bit(unsigned char *infile)
{
	if (bits_to_go == 0) {			/* Read the next byte if no	*/

		buffer2 = infile[nextchar];
		nextchar++;
		
		bits_to_go = 8;
	}
	/*
	 * Return the next bit
	 */
	bits_to_go -= 1;
	return((buffer2>>bits_to_go) & 1);
}

/*  ############################################################################  */
/* INPUT N BITS (N must be <= 8) */

static int input_nbits(unsigned char *infile, int n)
{
    /* AND mask for retreiving the right-most n bits */
    static int mask[9] = {0, 1, 3, 7, 15, 31, 63, 127, 255};

	if (bits_to_go < n) {
		/*
		 * need another byte's worth of bits
		 */

		buffer2 = (buffer2<<8) | (int) infile[nextchar];
		nextchar++;
		bits_to_go += 8;
	}
	/*
	 * now pick off the first n bits
	 */
	bits_to_go -= n;

        /* there was a slight gain in speed by replacing the following line */
/*	return( (buffer2>>bits_to_go) & ((1<<n)-1) ); */
	return( (buffer2>>bits_to_go) & (*(mask+n)) ); 
}
/*  ############################################################################  */
/* INPUT 4 BITS  */

static int input_nybble(unsigned char *infile)
{
	if (bits_to_go < 4) {
		/*
		 * need another byte's worth of bits
		 */

		buffer2 = (buffer2<<8) | (int) infile[nextchar];
		nextchar++;
		bits_to_go += 8;
	}
	/*
	 * now pick off the first 4 bits
	 */
	bits_to_go -= 4;

	return( (buffer2>>bits_to_go) & 15 ); 
}
/*  ############################################################################  */
/* INPUT array of 4 BITS  */

static int input_nnybble(unsigned char *infile, int n, unsigned char array[])
{
	/* copy n 4-bit nybbles from infile to the lower 4 bits of array */

int ii, kk, shift1, shift2;

/*  forcing byte alignment doesn;t help, and even makes it go slightly slower
if (bits_to_go != 8) input_nbits(infile, bits_to_go);
*/
	if (n == 1) {
		array[0] = input_nybble(infile);
		return(0);
	}
	
	if (bits_to_go == 8) {
		/*
		   already have 2 full nybbles in buffer2, so 
		   backspace the infile array to reuse last char
		*/
		nextchar--;
		bits_to_go = 0;
	}
	
	/* bits_to_go now has a value in the range 0 - 7.  After adding  */
	/* another byte, bits_to_go effectively will be in range 8 - 15 */	

	shift1 = bits_to_go + 4;   /* shift1 will be in range 4 - 11 */
	shift2 = bits_to_go;	   /* shift2 will be in range 0 -  7 */
	kk = 0;

	/* special case */
	if (bits_to_go == 0) 
	{
	    for (ii = 0; ii < n/2; ii++) {
		/*
		 * refill the buffer with next byte
		 */
		buffer2 = (buffer2<<8) | (int) infile[nextchar];
		nextchar++;
		array[kk]     = (int) ((buffer2>>4) & 15);
		array[kk + 1] = (int) ((buffer2) & 15);    /* no shift required */
		kk += 2;
	    }
	}
	else
	{
	    for (ii = 0; ii < n/2; ii++) {
		/*
		 * refill the buffer with next byte
		 */
		buffer2 = (buffer2<<8) | (int) infile[nextchar];
		nextchar++;
		array[kk]     = (int) ((buffer2>>shift1) & 15);
		array[kk + 1] = (int) ((buffer2>>shift2) & 15);
		kk += 2;
	    }
	}


	if (ii * 2 != n) {  /* have to read last odd byte */
		array[n-1] = input_nybble(infile);
	}

	return( (buffer2>>bits_to_go) & 15 ); 
}