nom.tam.fits.compression.algorithm.hcompress.HCompress Maven / Gradle / Ivy
package nom.tam.fits.compression.algorithm.hcompress;
import java.nio.ByteBuffer;
import java.nio.LongBuffer;
/*
* #%L
* nom.tam FITS library
* %%
* Copyright (C) 1996 - 2024 nom-tam-fits
* %%
* This is free and unencumbered software released into the public domain.
*
* Anyone is free to copy, modify, publish, use, compile, sell, or
* distribute this software, either in source code form or as a compiled
* binary, for any purpose, commercial or non-commercial, and by any
* means.
*
* In jurisdictions that recognize copyright laws, the author or authors
* of this software dedicate any and all copyright interest in the
* software to the public domain. We make this dedication for the benefit
* of the public at large and to the detriment of our heirs and
* successors. We intend this dedication to be an overt act of
* relinquishment in perpetuity of all present and future rights to this
* software under copyright law.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR
* OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
* ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
* OTHER DEALINGS IN THE SOFTWARE.
* #L%
*/
/**
*
* (for internal use) A hierarchical data compression algoritm, used by the Hubble Data Archive and the STScI
* Digital Sky Survey.
*
*
* The original compression code was written by Richard White at the STScI and included (ported to c and adapted) in
* cfitsio by William Pence, NASA/GSFC. That code was then ported to java by R. van Nieuwenhoven. Later it was massively
* refactored to harmonize the different compression algorithms and reduce the duplicate code pieces without obscuring
* the algorithm itself as far as possible. The original site for the algorithm is
*
*
* See http://www.stsci.edu /software/hcompress.html
*
*
* @author Richard White
* @author William Pence
* @author Richard van Nieuwenhoven
*
* @see HDecompress
* @see HCompressorOption
*/
@SuppressWarnings("javadoc")
public class HCompress {
private static final int HTRANS_START_MASK = -2;
protected static final double ROUNDING_HALF = 0.5;
protected static final int BITS_OF_1_BYTE = 8;
protected static final int BITS_OF_1_NYBBLE = 4;
protected static final int BYTE_MASK = 0xff;
protected static final int NYBBLE_MASK = 0xF;
/**
* to be refactored to a good name.
*/
private static final int N3 = 3;
private static final int[] BITS_MASK = {0, 1, 3, 7, 15, 31, 63, 127, 255};
/*
* Huffman code values and number of bits in each code
*/
private static final int[] CODE = {0x3e, 0x00, 0x01, 0x08, 0x02, 0x09, 0x1a, 0x1b, 0x03, 0x1c, 0x0a, 0x1d, 0x0b, 0x1e,
0x3f, 0x0c};
private static final byte[] CODE_MAGIC = {(byte) 0xDD, (byte) 0x99};
private static final int[] NCODE = {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
*/
private int bitbuffer;
/** Number of bits free in buffer */
private int bitsToGo2;
private int bitsToGo3;
/** Bits buffered for output */
private int buffer2;
private int b2i(boolean b) {
return b ? 1 : 0;
}
private int bufcopy(byte[] a, int n, byte[] buffer, int b, long 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]] << bitsToGo3;
bitsToGo3 += NCODE[a[i]];
if (bitsToGo3 >= BITS_OF_1_BYTE) {
buffer[b] = (byte) (bitbuffer & BYTE_MASK);
b++;
/*
* return warning code if we fill buffer
*/
if (b >= bmax) {
return b;
}
bitbuffer >>= BITS_OF_1_BYTE;
bitsToGo3 -= BITS_OF_1_BYTE;
}
}
}
return b;
}
protected void compress(long[] aa, int ny, int nx, int scale, ByteBuffer output) {
/*
* compress the input image using the H-compress algorithm a - input image tiledImageOperation 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 tiledImageOperation 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
*/
/* H-transform */
htrans(aa, nx, ny);
LongBuffer a = LongBuffer.wrap(aa);
/* digitize */
digitize(a, 0, nx, ny, scale);
/* encode and write to output tiledImageOperation */
encode(output, a, nx, ny, scale);
}
private LongBuffer copy(LongBuffer a, int i) {
a.position(i);
return a.slice();
}
private void digitize(LongBuffer a, int aOffset, int nx, int ny, long scale) {
/*
* round to multiple of scale
*/
if (scale <= 1) {
return;
}
long d = (scale + 1L) / 2L - 1L;
for (int index = 0; index < a.limit(); index++) {
long current = a.get(index);
a.put(index, (current > 0 ? current + d : current - d) / scale);
}
}
/**
* encode pixels.
*
* @param compressedBytes compressed data
* @param pixels pixels to compress
* @param nx image width dimension
* @param ny image height dimension
* @param nbitplanes Number of bit planes in quadrants
*/
private void doEncode(ByteBuffer compressedBytes, LongBuffer pixels, int nx, int ny, byte[] nbitplanes) {
int nx2 = (nx + 1) / 2;
int ny2 = (ny + 1) / 2;
/*
* Initialize bit output
*/
startOutputtingBits();
/*
* write out the bit planes for each quadrant
*/
qtreeEncode(compressedBytes, copy(pixels, 0), ny, nx2, ny2, nbitplanes[0]);
qtreeEncode(compressedBytes, copy(pixels, ny2), ny, nx2, ny / 2, nbitplanes[1]);
qtreeEncode(compressedBytes, copy(pixels, ny * nx2), ny, nx / 2, ny2, nbitplanes[1]);
qtreeEncode(compressedBytes, copy(pixels, ny * nx2 + ny2), ny, nx / 2, ny / 2, nbitplanes[2]);
/*
* Add zero as an EOF symbol
*/
outputNybble(compressedBytes, 0);
doneOutputtingBits(compressedBytes);
}
private void doneOutputtingBits(ByteBuffer outfile) {
if (bitsToGo2 < BITS_OF_1_BYTE) {
/* putc(buffer2< 0) {
/*
* positive element, put zero at end of buffer
*/
signbits[nsign] <<= 1;
bitsToGo--;
} else if (a.get(i) < 0) {
/*
* negative element, shift in a one
*/
signbits[nsign] <<= 1;
signbits[nsign] |= 1;
bitsToGo--;
/*
* replace a by absolute value
*/
a.put(i, -a.get(i));
}
if (bitsToGo == 0) {
/*
* filled up this byte, go to the next one
*/
bitsToGo = BITS_OF_1_BYTE;
nsign++;
signbits[nsign] = 0;
}
}
if (bitsToGo != BITS_OF_1_BYTE) {
/*
* some bits in last element move bits in last byte to bottom and increment nsign
*/
signbits[nsign] <<= bitsToGo;
nsign++;
}
/*
* calculate number of bit planes for 3 quadrants quadrant 0=bottom left, 1=bottom right or top left, 2=top
* right,
*/
for (int q = 0; q < N3; q++) {
vmax[q] = 0;
}
/*
* get maximum absolute value in each quadrant
*/
int nx2 = (nx + 1) / 2;
int ny2 = (ny + 1) / 2;
int j = 0; /* column counter */
int k = 0; /* row counter */
for (int i = 0; i < nel; i++) {
int q = (j >= ny2 ? 1 : 0) + (k >= nx2 ? 1 : 0);
if (vmax[q] < a.get(i)) {
vmax[q] = a.get(i);
}
if (++j >= ny) {
j = 0;
k++;
}
}
/*
* now calculate number of bits for each quadrant
*/
/* this is a more efficient way to do this, */
for (int q = 0; q < N3; q++) {
nbitplanes[q] = 0;
while (vmax[q] > 0) {
vmax[q] = vmax[q] >> 1;
nbitplanes[q]++;
}
}
/*
* write nbitplanes
*/
compressedBytes.put(nbitplanes, 0, nbitplanes.length);
/*
* write coded tiledImageOperation
*/
doEncode(compressedBytes, a, nx, ny, nbitplanes);
/*
* write sign bits
*/
if (nsign > 0) {
compressedBytes.put(signbits, 0, nsign);
}
return (int) noutchar;
}
private int htrans(long[] a, int nx, int ny) {
/*
* log2n is log2 of max(nx,ny) rounded up to next power of 2
*/
int nmax = nx > ny ? nx : ny;
int log2n = log2n(nmax);
if (nmax > 1 << log2n) {
log2n++;
}
/*
* get temporary storage for shuffling elements
*/
long[] tmp = new long[(nmax + 1) / 2];
/*
* set up rounding and shifting masks
*/
long shift = 0;
long mask = HTRANS_START_MASK;
long mask2 = mask << 1;
long prnd = 1;
long prnd2 = prnd << 1;
long nrnd2 = prnd2 - 1;
/*
* do log2n reductions We're indexing a as a 2-D tiledImageOperation with dimensions (nx,ny).
*/
int nxtop = nx;
int nytop = ny;
for (int k = 0; k < log2n; k++) {
int oddx = nxtop % 2;
int oddy = nytop % 2;
int i = 0;
for (; i < nxtop - oddx; i += 2) {
int s00 = i * ny; /* s00 is index of a[i,j] */
int s10 = s00 + ny; /* s10 is index of a[i+1,j] */
for (int j = 0; j < nytop - oddy; j += 2) {
/*
* Divide h0,hx,hy,hc by 2 (1 the first time through).
*/
long h0 = a[s10 + 1] + a[s10] + a[s00 + 1] + a[s00] >> shift;
long hx = a[s10 + 1] + a[s10] - a[s00 + 1] - a[s00] >> shift;
long hy = a[s10 + 1] - a[s10] + a[s00 + 1] - a[s00] >> shift;
long 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 != 0) {
/*
* do last element in row if row length is odd s00+1, s10+1 are off edge
*/
long h0 = a[s10] + a[s00] << 1 - shift;
long 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++;
s10++;
}
}
if (oddx != 0) {
/*
* do last row if column length is odd s10, s10+1 are off edge
*/
int s00 = i * ny;
for (int j = 0; j < nytop - oddy; j += 2) {
long h0 = a[s00 + 1] + a[s00] << 1 - shift;
long 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 != 0) {
/*
* do corner element if both row and column lengths are odd s00+1, s10, s10+1 are off edge
*/
long h0 = a[s00] << 2 - shift;
a[s00] = (h0 >= 0 ? h0 + prnd2 : h0 + nrnd2) & mask2;
}
}
/*
* now shuffle in each dimension to group coefficients by order
*/
// achtung eigenlich pointer nach a
for (i = 0; i < nxtop; i++) {
shuffle(a, ny * i, nytop, 1, tmp);
}
for (int 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;
}
return 0;
}
private int log2n(int nqmax) {
return (int) (Math.log(nqmax) / Math.log(2.0) + ROUNDING_HALF);
}
private void outputNbits(ByteBuffer outfile, int bits, int n) {
/* AND mask for the right-most n bits */
/*
* insert bits at end of buffer
*/
buffer2 <<= n;
/* buffer2 |= ( bits & ((1<> -bitsToGo2 & BYTE_MASK));
bitsToGo2 += BITS_OF_1_BYTE;
}
}
private void outputNnybble(ByteBuffer outfile, int n, byte[] array) {
/*
* pack the 4 lower bits in each element of the tiledImageOperation into the outfile tiledImageOperation
*/
int ii, jj, kk = 0, shift;
if (n == 1) {
outputNybble(outfile, 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 (bitsToGo2 <= BITS_OF_1_NYBBLE) {
/* just room for 1 nybble; write it out separately */
outputNybble(outfile, array[0]);
kk++; /* index to next tiledImageOperation element */
if (n == 2) {
// only 1 more nybble to write out
outputNybble(outfile, array[1]);
return;
}
}
/* bits_to_go2 is now in the range 5 - 8 */
shift = BITS_OF_1_BYTE - bitsToGo2;
/*
* now write out pairs of nybbles; this does not affect value of bits_to_go2
*/
jj = (n - kk) / 2;
if (bitsToGo2 == BITS_OF_1_BYTE) {
/* special case if nybbles are aligned on byte boundary */
/* this actually seems to make very little difference in speed */
buffer2 = 0;
for (ii = 0; ii < jj; ii++) {
outfile.put((byte) ((array[kk] & NYBBLE_MASK) << BITS_OF_1_NYBBLE | array[kk + 1] & NYBBLE_MASK));
kk += 2;
}
} else {
for (ii = 0; ii < jj; ii++) {
buffer2 = buffer2 << BITS_OF_1_BYTE | (array[kk] & NYBBLE_MASK) << BITS_OF_1_NYBBLE
| array[kk + 1] & NYBBLE_MASK;
kk += 2;
/*
* buffer2 full, put out top 8 bits
*/
outfile.put((byte) (buffer2 >> shift & BYTE_MASK));
}
}
/* write out last odd nybble, if present */
if (kk != n) {
outputNybble(outfile, array[n - 1]);
}
}
private void outputNybble(ByteBuffer outfile, int bits) {
/*
* insert 4 bits at end of buffer
*/
buffer2 = buffer2 << BITS_OF_1_NYBBLE | bits & NYBBLE_MASK;
bitsToGo2 -= BITS_OF_1_NYBBLE;
if (bitsToGo2 <= 0) {
/*
* buffer2 full, put out top 8 bits
*/
outfile.put((byte) (buffer2 >> -bitsToGo2 & BYTE_MASK));
bitsToGo2 += BITS_OF_1_BYTE;
}
}
/**
* macros to write out 4-bit nybble, Huffman code for this value
*/
private int qtreeEncode(ByteBuffer outfile, LongBuffer 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, nqmax, nqx2, nqy2, nx, ny;
long bmax;
byte[] scratch, buffer;
/*
* log2n is log2 of max(nqx,nqy) rounded up to next power of 2
*/
nqmax = nqx > nqy ? nqx : nqy;
log2n = log2n(nqmax);
if (nqmax > 1 << log2n) {
log2n++;
}
/*
* 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 tiledImageOperation 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 = new byte[(int) (2 * bmax)];
buffer = new byte[(int) bmax];
/*
* now encode each bit plane, starting with the top
*/
bitplane_done: for (bit = nbitplanes - 1; bit >= 0; bit--) {
/*
* initial bit buffer
*/
b = 0;
bitbuffer = 0;
bitsToGo3 = 0;
/*
* on first pass copy A to scratch tiledImageOperation
*/
qtreeOnebit(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
*/
b = bufcopy(scratch, nx * ny, buffer, b, bmax);
if (b >= bmax) {
/*
* quadtree is expanding data, change warning code and just fill buffer with bit-map
*/
writeBdirect(outfile, a, n, nqx, nqy, scratch, bit);
continue bitplane_done;
}
/*
* do log2n reductions
*/
for (k = 1; k < log2n; k++) {
qtreeReduce(scratch, ny, nx, ny, scratch);
nx = nx + 1 >> 1;
ny = ny + 1 >> 1;
b = bufcopy(scratch, nx * ny, buffer, b, bmax);
if (b >= bmax) {
writeBdirect(outfile, a, n, nqx, nqy, scratch, bit);
continue bitplane_done;
}
}
/*
* OK, we've got the code in buffer Write quadtree warning code, then write buffer in reverse order
*/
outputNybble(outfile, NYBBLE_MASK);
if (b == 0) {
if (bitsToGo3 > 0) {
/*
* put out the last few bits
*/
outputNbits(outfile, bitbuffer & (1 << bitsToGo3) - 1, bitsToGo3);
} else {
/*
* have to write a zero nybble if there are no 1's in tiledImageOperation
*/
outputNbits(outfile, CODE[0], NCODE[0]);
}
} else {
if (bitsToGo3 > 0) {
/*
* put out the last few bits
*/
outputNbits(outfile, bitbuffer & (1 << bitsToGo3) - 1, bitsToGo3);
}
for (i = b - 1; i >= 0; i--) {
outputNbits(outfile, buffer[i], BITS_OF_1_BYTE);
}
}
}
return 0;
}
private void qtreeOnebit(LongBuffer a, int n, int nx, int ny, byte[] b, int bit) {
int i, j, k;
long b0, b1, b2, b3;
int s10, s00;
/*
* use selected bit to get amount to shift
*/
b0 = 1L << bit;
b1 = b0 << 1;
b2 = b1 << 1;
b3 = b2 << 1;
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) {
b[k] = (byte) ((a.get(s10 + 1) & b0 //
| a.get(s10) << 1 & b1 //
| a.get(s00 + 1) << 2 & b2 //
| a.get(s00) << N3 & b3) >> bit);
k++;
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] = (byte) ((a.get(s10) << 1 & b1 | a.get(s00) << N3 & b3) >> bit);
k++;
}
}
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] = (byte) ((a.get(s00 + 1) << 2 & b2 | a.get(s00) << N3 & b3) >> bit);
k++;
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] = (byte) ((a.get(s00) << N3 & b3) >> bit);
k++;
}
}
}
private void qtreeReduce(byte[] a, int n, int nx, int ny, byte[] 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] = (byte) (b2i(a[s10 + 1] != 0) | b2i(a[s10] != 0) << 1 | b2i(a[s00 + 1] != 0) << 2
| b2i(a[s00] != 0) << N3);
k++;
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] = (byte) (b2i(a[s10] != 0) << 1 | b2i(a[s00] != 0) << N3);
k++;
}
}
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] = (byte) (b2i(a[s00 + 1] != 0) << 2 | b2i(a[s00] != 0) << N3);
k++;
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] = (byte) (b2i(a[s00] != 0) << N3);
k++;
}
}
}
private void shuffle(long[] a, int aOffset, int n, int n2, long[] tmp) {
/*
* int a[]; tiledImageOperation to shuffle int n; number of elements to shuffle int n2; second dimension int
* tmp[]; scratch storage
*/
int i;
long[] p1, p2, pt;
int p1Offset;
int ptOffset;
int p2Offset;
/*
* copy odd elements to tmp
*/
pt = tmp;
ptOffset = 0;
p1 = a;
p1Offset = aOffset + n2;
for (i = 1; i < n; i += 2) {
pt[ptOffset] = p1[p1Offset];
ptOffset++;
p1Offset += n2 + n2;
}
/*
* compress even elements into first half of A
*/
p1 = a;
p1Offset = aOffset + n2;
p2 = a;
p2Offset = aOffset + n2 + n2;
for (i = 2; i < n; i += 2) {
p1[p1Offset] = p2[p2Offset];
p1Offset += n2;
p2Offset += n2 + n2;
}
/*
* put odd elements into 2nd half
*/
pt = tmp;
ptOffset = 0;
for (i = 1; i < n; i += 2) {
p1[p1Offset] = pt[ptOffset];
p1Offset += n2;
ptOffset++;
}
}
private void startOutputtingBits() {
buffer2 = 0; /* Buffer is empty to start */
bitsToGo2 = BITS_OF_1_BYTE; /* with */
}
private void writeBdirect(ByteBuffer outfile, LongBuffer a, int n, int nqx, int nqy, byte[] scratch, int bit) {
/*
* Write the direct bitmap warning code
*/
outputNybble(outfile, 0x0);
/*
* Copy A to scratch tiledImageOperation (again!), packing 4 bits/nybble
*/
qtreeOnebit(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]); }
*/
outputNnybble(outfile, (nqx + 1) / 2 * ((nqy + 1) / 2), scratch);
}
}