1 /* #########################################################################
2 These routines to apply the H-compress decompression algorithm to a 2-D Fits
3 image were written by R. White at the STScI and were obtained from the STScI at
4 http://www.stsci.edu/software/hcompress.html
5
6 This source file is a concatination of the following sources files in the
7 original distribution
8 hinv.c
9 hsmooth.c
10 undigitize.c
11 decode.c
12 dodecode.c
13 qtree_decode.c
14 qread.c
15 bit_input.c
16
17
18 The following modifications have been made to the original code:
19
20 - commented out redundant "include" statements
21 - added the nextchar global variable
22 - changed all the 'extern' declarations to 'static', since all the routines are in
23 the same source file
24 - changed the first parameter in decode (and in lower level routines from a file stream
25 to a char array
26 - modified the myread routine, and lower level byte reading routines, to copy
27 the input bytes to a char array, instead of reading them from a file stream
28 - changed the function declarations to the more modern ANSI C style
29 - changed calls to printf and perror to call the CFITSIO ffpmsg routine
30 - replace "exit" statements with "return" statements
31
32 ############################################################################ */
33
34 #include <stdio.h>
35 #include <math.h>
36 #include <stdlib.h>
37 #include <string.h>
38 #include "fitsio2.h"
39
40 /* WDP added test to see if min and max are already defined */
41 #ifndef min
42 #define min(a,b) (((a)<(b))?(a):(b))
43 #endif
44 #ifndef max
45 #define max(a,b) (((a)>(b))?(a):(b))
46 #endif
47
48 static long nextchar;
49
50 static int decode(unsigned char *infile, int *a, int *nx, int *ny, int *scale);
51 static int decode64(unsigned char *infile, LONGLONG *a, int *nx, int *ny, int *scale);
52 static int hinv(int a[], int nx, int ny, int smooth ,int scale);
53 static int hinv64(LONGLONG a[], int nx, int ny, int smooth ,int scale);
54 static void undigitize(int a[], int nx, int ny, int scale);
55 static void undigitize64(LONGLONG a[], int nx, int ny, int scale);
56 static void unshuffle(int a[], int n, int n2, int tmp[]);
57 static void unshuffle64(LONGLONG a[], int n, int n2, LONGLONG tmp[]);
58 static void hsmooth(int a[], int nxtop, int nytop, int ny, int scale);
59 static void hsmooth64(LONGLONG a[], int nxtop, int nytop, int ny, int scale);
60 static void qread(unsigned char *infile,char *a, int n);
61 static int readint(unsigned char *infile);
62 static LONGLONG readlonglong(unsigned char *infile);
63 static int dodecode(unsigned char *infile, int a[], int nx, int ny, unsigned char nbitplanes[3]);
64 static int dodecode64(unsigned char *infile, LONGLONG a[], int nx, int ny, unsigned char nbitplanes[3]);
65 static int qtree_decode(unsigned char *infile, int a[], int n, int nqx, int nqy, int nbitplanes);
66 static int qtree_decode64(unsigned char *infile, LONGLONG a[], int n, int nqx, int nqy, int nbitplanes);
67 static void start_inputing_bits(void);
68 static int input_bit(unsigned char *infile);
69 static int input_nbits(unsigned char *infile, int n);
70 /* make input_nybble a separate routine, for added effiency */
71 /* #define input_nybble(infile) input_nbits(infile,4) */
72 static int input_nybble(unsigned char *infile);
73 static int input_nnybble(unsigned char *infile, int n, unsigned char *array);
74
75 static void qtree_expand(unsigned char *infile, unsigned char a[], int nx, int ny, unsigned char b[]);
76 static void qtree_bitins(unsigned char a[], int nx, int ny, int b[], int n, int bit);
77 static void qtree_bitins64(unsigned char a[], int nx, int ny, LONGLONG b[], int n, int bit);
78 static void qtree_copy(unsigned char a[], int nx, int ny, unsigned char b[], int n);
79 static void read_bdirect(unsigned char *infile, int a[], int n, int nqx, int nqy, unsigned char scratch[], int bit);
80 static void read_bdirect64(unsigned char *infile, LONGLONG a[], int n, int nqx, int nqy, unsigned char scratch[], int bit);
81 static int input_huffman(unsigned char *infile);
82
83 /* ---------------------------------------------------------------------- */
fits_hdecompress(unsigned char * input,int smooth,int * a,int * ny,int * nx,int * scale,int * status)84 int fits_hdecompress(unsigned char *input, int smooth, int *a, int *ny, int *nx,
85 int *scale, int *status)
86 {
87 /*
88 decompress the input byte stream using the H-compress algorithm
89
90 input - input array of compressed bytes
91 a - pre-allocated array to hold the output uncompressed image
92 nx - returned X axis size
93 ny - returned Y axis size
94
95 NOTE: the nx and ny dimensions as defined within this code are reversed from
96 the usual FITS notation. ny is the fastest varying dimension, which is
97 usually considered the X axis in the FITS image display
98
99 */
100 int stat;
101
102 if (*status > 0) return(*status);
103
104 /* decode the input array */
105
106 FFLOCK; /* decode uses the nextchar global variable */
107 stat = decode(input, a, nx, ny, scale);
108 FFUNLOCK;
109
110 *status = stat;
111 if (stat) return(*status);
112
113 /*
114 * Un-Digitize
115 */
116 undigitize(a, *nx, *ny, *scale);
117
118 /*
119 * Inverse H-transform
120 */
121 stat = hinv(a, *nx, *ny, smooth, *scale);
122 *status = stat;
123
124 return(*status);
125 }
126 /* ---------------------------------------------------------------------- */
fits_hdecompress64(unsigned char * input,int smooth,LONGLONG * a,int * ny,int * nx,int * scale,int * status)127 int fits_hdecompress64(unsigned char *input, int smooth, LONGLONG *a, int *ny, int *nx,
128 int *scale, int *status)
129 {
130 /*
131 decompress the input byte stream using the H-compress algorithm
132
133 input - input array of compressed bytes
134 a - pre-allocated array to hold the output uncompressed image
135 nx - returned X axis size
136 ny - returned Y axis size
137
138 NOTE: the nx and ny dimensions as defined within this code are reversed from
139 the usual FITS notation. ny is the fastest varying dimension, which is
140 usually considered the X axis in the FITS image display
141
142 */
143 int stat, *iarray, ii, nval;
144
145 if (*status > 0) return(*status);
146
147 /* decode the input array */
148
149 FFLOCK; /* decode uses the nextchar global variable */
150 stat = decode64(input, a, nx, ny, scale);
151 FFUNLOCK;
152
153 *status = stat;
154 if (stat) return(*status);
155
156 /*
157 * Un-Digitize
158 */
159 undigitize64(a, *nx, *ny, *scale);
160
161 /*
162 * Inverse H-transform
163 */
164 stat = hinv64(a, *nx, *ny, smooth, *scale);
165
166 *status = stat;
167
168 /* pack the I*8 values back into an I*4 array */
169 iarray = (int *) a;
170 nval = (*nx) * (*ny);
171
172 for (ii = 0; ii < nval; ii++)
173 iarray[ii] = (int) a[ii];
174
175 return(*status);
176 }
177
178 /* ############################################################################ */
179 /* ############################################################################ */
180
181 /* Copyright (c) 1993 Association of Universities for Research
182 * in Astronomy. All rights reserved. Produced under National
183 * Aeronautics and Space Administration Contract No. NAS5-26555.
184 */
185 /* hinv.c Inverse H-transform of NX x NY integer image
186 *
187 * Programmer: R. White Date: 23 July 1993
188 */
189
190 /* ############################################################################ */
191 static int
hinv(int a[],int nx,int ny,int smooth,int scale)192 hinv(int a[], int nx, int ny, int smooth ,int scale)
193 /*
194 int smooth; 0 for no smoothing, else smooth during inversion
195 int scale; used if smoothing is specified
196 */
197 {
198 int nmax, log2n, i, j, k;
199 int nxtop,nytop,nxf,nyf,c;
200 int oddx,oddy;
201 int shift, bit0, bit1, bit2, mask0, mask1, mask2,
202 prnd0, prnd1, prnd2, nrnd0, nrnd1, nrnd2, lowbit0, lowbit1;
203 int h0, hx, hy, hc;
204 int s10, s00;
205 int *tmp;
206
207 /*
208 * log2n is log2 of max(nx,ny) rounded up to next power of 2
209 */
210 nmax = (nx>ny) ? nx : ny;
211 log2n = (int) (log((float) nmax)/log(2.0)+0.5);
212 if ( nmax > (1<<log2n) ) {
213 log2n += 1;
214 }
215 /*
216 * get temporary storage for shuffling elements
217 */
218 tmp = (int *) malloc(((nmax+1)/2)*sizeof(int));
219 if (tmp == (int *) NULL) {
220 ffpmsg("hinv: insufficient memory");
221 return(DATA_DECOMPRESSION_ERR);
222 }
223 /*
224 * set up masks, rounding parameters
225 */
226 shift = 1;
227 bit0 = 1 << (log2n - 1);
228 bit1 = bit0 << 1;
229 bit2 = bit0 << 2;
230 mask0 = -bit0;
231 mask1 = mask0 << 1;
232 mask2 = mask0 << 2;
233 prnd0 = bit0 >> 1;
234 prnd1 = bit1 >> 1;
235 prnd2 = bit2 >> 1;
236 nrnd0 = prnd0 - 1;
237 nrnd1 = prnd1 - 1;
238 nrnd2 = prnd2 - 1;
239 /*
240 * round h0 to multiple of bit2
241 */
242 a[0] = (a[0] + ((a[0] >= 0) ? prnd2 : nrnd2)) & mask2;
243 /*
244 * do log2n expansions
245 *
246 * We're indexing a as a 2-D array with dimensions (nx,ny).
247 */
248 nxtop = 1;
249 nytop = 1;
250 nxf = nx;
251 nyf = ny;
252 c = 1<<log2n;
253 for (k = log2n-1; k>=0; k--) {
254 /*
255 * this somewhat cryptic code generates the sequence
256 * ntop[k-1] = (ntop[k]+1)/2, where ntop[log2n] = n
257 */
258 c = c>>1;
259 nxtop = nxtop<<1;
260 nytop = nytop<<1;
261 if (nxf <= c) { nxtop -= 1; } else { nxf -= c; }
262 if (nyf <= c) { nytop -= 1; } else { nyf -= c; }
263 /*
264 * double shift and fix nrnd0 (because prnd0=0) on last pass
265 */
266 if (k == 0) {
267 nrnd0 = 0;
268 shift = 2;
269 }
270 /*
271 * unshuffle in each dimension to interleave coefficients
272 */
273 for (i = 0; i<nxtop; i++) {
274 unshuffle(&a[ny*i],nytop,1,tmp);
275 }
276 for (j = 0; j<nytop; j++) {
277 unshuffle(&a[j],nxtop,ny,tmp);
278 }
279 /*
280 * smooth by interpolating coefficients if SMOOTH != 0
281 */
282 if (smooth) hsmooth(a,nxtop,nytop,ny,scale);
283 oddx = nxtop % 2;
284 oddy = nytop % 2;
285 for (i = 0; i<nxtop-oddx; i += 2) {
286 s00 = ny*i; /* s00 is index of a[i,j] */
287 s10 = s00+ny; /* s10 is index of a[i+1,j] */
288 for (j = 0; j<nytop-oddy; j += 2) {
289 h0 = a[s00 ];
290 hx = a[s10 ];
291 hy = a[s00+1];
292 hc = a[s10+1];
293 /*
294 * round hx and hy to multiple of bit1, hc to multiple of bit0
295 * h0 is already a multiple of bit2
296 */
297 hx = (hx + ((hx >= 0) ? prnd1 : nrnd1)) & mask1;
298 hy = (hy + ((hy >= 0) ? prnd1 : nrnd1)) & mask1;
299 hc = (hc + ((hc >= 0) ? prnd0 : nrnd0)) & mask0;
300 /*
301 * propagate bit0 of hc to hx,hy
302 */
303 lowbit0 = hc & bit0;
304 hx = (hx >= 0) ? (hx - lowbit0) : (hx + lowbit0);
305 hy = (hy >= 0) ? (hy - lowbit0) : (hy + lowbit0);
306 /*
307 * Propagate bits 0 and 1 of hc,hx,hy to h0.
308 * This could be simplified if we assume h0>0, but then
309 * the inversion would not be lossless for images with
310 * negative pixels.
311 */
312 lowbit1 = (hc ^ hx ^ hy) & bit1;
313 h0 = (h0 >= 0)
314 ? (h0 + lowbit0 - lowbit1)
315 : (h0 + ((lowbit0 == 0) ? lowbit1 : (lowbit0-lowbit1)));
316 /*
317 * Divide sums by 2 (4 last time)
318 */
319 a[s10+1] = (h0 + hx + hy + hc) >> shift;
320 a[s10 ] = (h0 + hx - hy - hc) >> shift;
321 a[s00+1] = (h0 - hx + hy - hc) >> shift;
322 a[s00 ] = (h0 - hx - hy + hc) >> shift;
323 s00 += 2;
324 s10 += 2;
325 }
326 if (oddy) {
327 /*
328 * do last element in row if row length is odd
329 * s00+1, s10+1 are off edge
330 */
331 h0 = a[s00 ];
332 hx = a[s10 ];
333 hx = ((hx >= 0) ? (hx+prnd1) : (hx+nrnd1)) & mask1;
334 lowbit1 = hx & bit1;
335 h0 = (h0 >= 0) ? (h0 - lowbit1) : (h0 + lowbit1);
336 a[s10 ] = (h0 + hx) >> shift;
337 a[s00 ] = (h0 - hx) >> shift;
338 }
339 }
340 if (oddx) {
341 /*
342 * do last row if column length is odd
343 * s10, s10+1 are off edge
344 */
345 s00 = ny*i;
346 for (j = 0; j<nytop-oddy; j += 2) {
347 h0 = a[s00 ];
348 hy = a[s00+1];
349 hy = ((hy >= 0) ? (hy+prnd1) : (hy+nrnd1)) & mask1;
350 lowbit1 = hy & bit1;
351 h0 = (h0 >= 0) ? (h0 - lowbit1) : (h0 + lowbit1);
352 a[s00+1] = (h0 + hy) >> shift;
353 a[s00 ] = (h0 - hy) >> shift;
354 s00 += 2;
355 }
356 if (oddy) {
357 /*
358 * do corner element if both row and column lengths are odd
359 * s00+1, s10, s10+1 are off edge
360 */
361 h0 = a[s00 ];
362 a[s00 ] = h0 >> shift;
363 }
364 }
365 /*
366 * divide all the masks and rounding values by 2
367 */
368 bit2 = bit1;
369 bit1 = bit0;
370 bit0 = bit0 >> 1;
371 mask1 = mask0;
372 mask0 = mask0 >> 1;
373 prnd1 = prnd0;
374 prnd0 = prnd0 >> 1;
375 nrnd1 = nrnd0;
376 nrnd0 = prnd0 - 1;
377 }
378 free(tmp);
379 return(0);
380 }
381 /* ############################################################################ */
382 static int
hinv64(LONGLONG a[],int nx,int ny,int smooth,int scale)383 hinv64(LONGLONG a[], int nx, int ny, int smooth ,int scale)
384 /*
385 int smooth; 0 for no smoothing, else smooth during inversion
386 int scale; used if smoothing is specified
387 */
388 {
389 int nmax, log2n, i, j, k;
390 int nxtop,nytop,nxf,nyf,c;
391 int oddx,oddy;
392 int shift;
393 LONGLONG mask0, mask1, mask2, prnd0, prnd1, prnd2, bit0, bit1, bit2;
394 LONGLONG nrnd0, nrnd1, nrnd2, lowbit0, lowbit1;
395 LONGLONG h0, hx, hy, hc;
396 int s10, s00;
397 LONGLONG *tmp;
398
399 /*
400 * log2n is log2 of max(nx,ny) rounded up to next power of 2
401 */
402 nmax = (nx>ny) ? nx : ny;
403 log2n = (int) (log((float) nmax)/log(2.0)+0.5);
404 if ( nmax > (1<<log2n) ) {
405 log2n += 1;
406 }
407 /*
408 * get temporary storage for shuffling elements
409 */
410 tmp = (LONGLONG *) malloc(((nmax+1)/2)*sizeof(LONGLONG));
411 if (tmp == (LONGLONG *) NULL) {
412 ffpmsg("hinv64: insufficient memory");
413 return(DATA_DECOMPRESSION_ERR);
414 }
415 /*
416 * set up masks, rounding parameters
417 */
418 shift = 1;
419 bit0 = ((LONGLONG) 1) << (log2n - 1);
420 bit1 = bit0 << 1;
421 bit2 = bit0 << 2;
422 mask0 = -bit0;
423 mask1 = mask0 << 1;
424 mask2 = mask0 << 2;
425 prnd0 = bit0 >> 1;
426 prnd1 = bit1 >> 1;
427 prnd2 = bit2 >> 1;
428 nrnd0 = prnd0 - 1;
429 nrnd1 = prnd1 - 1;
430 nrnd2 = prnd2 - 1;
431 /*
432 * round h0 to multiple of bit2
433 */
434 a[0] = (a[0] + ((a[0] >= 0) ? prnd2 : nrnd2)) & mask2;
435 /*
436 * do log2n expansions
437 *
438 * We're indexing a as a 2-D array with dimensions (nx,ny).
439 */
440 nxtop = 1;
441 nytop = 1;
442 nxf = nx;
443 nyf = ny;
444 c = 1<<log2n;
445 for (k = log2n-1; k>=0; k--) {
446 /*
447 * this somewhat cryptic code generates the sequence
448 * ntop[k-1] = (ntop[k]+1)/2, where ntop[log2n] = n
449 */
450 c = c>>1;
451 nxtop = nxtop<<1;
452 nytop = nytop<<1;
453 if (nxf <= c) { nxtop -= 1; } else { nxf -= c; }
454 if (nyf <= c) { nytop -= 1; } else { nyf -= c; }
455 /*
456 * double shift and fix nrnd0 (because prnd0=0) on last pass
457 */
458 if (k == 0) {
459 nrnd0 = 0;
460 shift = 2;
461 }
462 /*
463 * unshuffle in each dimension to interleave coefficients
464 */
465 for (i = 0; i<nxtop; i++) {
466 unshuffle64(&a[ny*i],nytop,1,tmp);
467 }
468 for (j = 0; j<nytop; j++) {
469 unshuffle64(&a[j],nxtop,ny,tmp);
470 }
471 /*
472 * smooth by interpolating coefficients if SMOOTH != 0
473 */
474 if (smooth) hsmooth64(a,nxtop,nytop,ny,scale);
475 oddx = nxtop % 2;
476 oddy = nytop % 2;
477 for (i = 0; i<nxtop-oddx; i += 2) {
478 s00 = ny*i; /* s00 is index of a[i,j] */
479 s10 = s00+ny; /* s10 is index of a[i+1,j] */
480 for (j = 0; j<nytop-oddy; j += 2) {
481 h0 = a[s00 ];
482 hx = a[s10 ];
483 hy = a[s00+1];
484 hc = a[s10+1];
485 /*
486 * round hx and hy to multiple of bit1, hc to multiple of bit0
487 * h0 is already a multiple of bit2
488 */
489 hx = (hx + ((hx >= 0) ? prnd1 : nrnd1)) & mask1;
490 hy = (hy + ((hy >= 0) ? prnd1 : nrnd1)) & mask1;
491 hc = (hc + ((hc >= 0) ? prnd0 : nrnd0)) & mask0;
492 /*
493 * propagate bit0 of hc to hx,hy
494 */
495 lowbit0 = hc & bit0;
496 hx = (hx >= 0) ? (hx - lowbit0) : (hx + lowbit0);
497 hy = (hy >= 0) ? (hy - lowbit0) : (hy + lowbit0);
498 /*
499 * Propagate bits 0 and 1 of hc,hx,hy to h0.
500 * This could be simplified if we assume h0>0, but then
501 * the inversion would not be lossless for images with
502 * negative pixels.
503 */
504 lowbit1 = (hc ^ hx ^ hy) & bit1;
505 h0 = (h0 >= 0)
506 ? (h0 + lowbit0 - lowbit1)
507 : (h0 + ((lowbit0 == 0) ? lowbit1 : (lowbit0-lowbit1)));
508 /*
509 * Divide sums by 2 (4 last time)
510 */
511 a[s10+1] = (h0 + hx + hy + hc) >> shift;
512 a[s10 ] = (h0 + hx - hy - hc) >> shift;
513 a[s00+1] = (h0 - hx + hy - hc) >> shift;
514 a[s00 ] = (h0 - hx - hy + hc) >> shift;
515 s00 += 2;
516 s10 += 2;
517 }
518 if (oddy) {
519 /*
520 * do last element in row if row length is odd
521 * s00+1, s10+1 are off edge
522 */
523 h0 = a[s00 ];
524 hx = a[s10 ];
525 hx = ((hx >= 0) ? (hx+prnd1) : (hx+nrnd1)) & mask1;
526 lowbit1 = hx & bit1;
527 h0 = (h0 >= 0) ? (h0 - lowbit1) : (h0 + lowbit1);
528 a[s10 ] = (h0 + hx) >> shift;
529 a[s00 ] = (h0 - hx) >> shift;
530 }
531 }
532 if (oddx) {
533 /*
534 * do last row if column length is odd
535 * s10, s10+1 are off edge
536 */
537 s00 = ny*i;
538 for (j = 0; j<nytop-oddy; j += 2) {
539 h0 = a[s00 ];
540 hy = a[s00+1];
541 hy = ((hy >= 0) ? (hy+prnd1) : (hy+nrnd1)) & mask1;
542 lowbit1 = hy & bit1;
543 h0 = (h0 >= 0) ? (h0 - lowbit1) : (h0 + lowbit1);
544 a[s00+1] = (h0 + hy) >> shift;
545 a[s00 ] = (h0 - hy) >> shift;
546 s00 += 2;
547 }
548 if (oddy) {
549 /*
550 * do corner element if both row and column lengths are odd
551 * s00+1, s10, s10+1 are off edge
552 */
553 h0 = a[s00 ];
554 a[s00 ] = h0 >> shift;
555 }
556 }
557 /*
558 * divide all the masks and rounding values by 2
559 */
560 bit2 = bit1;
561 bit1 = bit0;
562 bit0 = bit0 >> 1;
563 mask1 = mask0;
564 mask0 = mask0 >> 1;
565 prnd1 = prnd0;
566 prnd0 = prnd0 >> 1;
567 nrnd1 = nrnd0;
568 nrnd0 = prnd0 - 1;
569 }
570 free(tmp);
571 return(0);
572 }
573
574 /* ############################################################################ */
575 static void
unshuffle(int a[],int n,int n2,int tmp[])576 unshuffle(int a[], int n, int n2, int tmp[])
577 /*
578 int a[]; array to shuffle
579 int n; number of elements to shuffle
580 int n2; second dimension
581 int tmp[]; scratch storage
582 */
583 {
584 int i;
585 int nhalf;
586 int *p1, *p2, *pt;
587
588 /*
589 * copy 2nd half of array to tmp
590 */
591 nhalf = (n+1)>>1;
592 pt = tmp;
593 p1 = &a[n2*nhalf]; /* pointer to a[i] */
594 for (i=nhalf; i<n; i++) {
595 *pt = *p1;
596 p1 += n2;
597 pt += 1;
598 }
599 /*
600 * distribute 1st half of array to even elements
601 */
602 p2 = &a[ n2*(nhalf-1) ]; /* pointer to a[i] */
603 p1 = &a[(n2*(nhalf-1))<<1]; /* pointer to a[2*i] */
604 for (i=nhalf-1; i >= 0; i--) {
605 *p1 = *p2;
606 p2 -= n2;
607 p1 -= (n2+n2);
608 }
609 /*
610 * now distribute 2nd half of array (in tmp) to odd elements
611 */
612 pt = tmp;
613 p1 = &a[n2]; /* pointer to a[i] */
614 for (i=1; i<n; i += 2) {
615 *p1 = *pt;
616 p1 += (n2+n2);
617 pt += 1;
618 }
619 }
620 /* ############################################################################ */
621 static void
unshuffle64(LONGLONG a[],int n,int n2,LONGLONG tmp[])622 unshuffle64(LONGLONG a[], int n, int n2, LONGLONG tmp[])
623 /*
624 LONGLONG a[]; array to shuffle
625 int n; number of elements to shuffle
626 int n2; second dimension
627 LONGLONG tmp[]; scratch storage
628 */
629 {
630 int i;
631 int nhalf;
632 LONGLONG *p1, *p2, *pt;
633
634 /*
635 * copy 2nd half of array to tmp
636 */
637 nhalf = (n+1)>>1;
638 pt = tmp;
639 p1 = &a[n2*nhalf]; /* pointer to a[i] */
640 for (i=nhalf; i<n; i++) {
641 *pt = *p1;
642 p1 += n2;
643 pt += 1;
644 }
645 /*
646 * distribute 1st half of array to even elements
647 */
648 p2 = &a[ n2*(nhalf-1) ]; /* pointer to a[i] */
649 p1 = &a[(n2*(nhalf-1))<<1]; /* pointer to a[2*i] */
650 for (i=nhalf-1; i >= 0; i--) {
651 *p1 = *p2;
652 p2 -= n2;
653 p1 -= (n2+n2);
654 }
655 /*
656 * now distribute 2nd half of array (in tmp) to odd elements
657 */
658 pt = tmp;
659 p1 = &a[n2]; /* pointer to a[i] */
660 for (i=1; i<n; i += 2) {
661 *p1 = *pt;
662 p1 += (n2+n2);
663 pt += 1;
664 }
665 }
666
667 /* ############################################################################ */
668 /* ############################################################################ */
669
670 /* Copyright (c) 1993 Association of Universities for Research
671 * in Astronomy. All rights reserved. Produced under National
672 * Aeronautics and Space Administration Contract No. NAS5-26555.
673 */
674 /* hsmooth.c Smooth H-transform image by adjusting coefficients toward
675 * interpolated values
676 *
677 * Programmer: R. White Date: 13 April 1992
678 */
679
680 /* ############################################################################ */
681 static void
hsmooth(int a[],int nxtop,int nytop,int ny,int scale)682 hsmooth(int a[], int nxtop, int nytop, int ny, int scale)
683 /*
684 int a[]; array of H-transform coefficients
685 int nxtop,nytop; size of coefficient block to use
686 int ny; actual 1st dimension of array
687 int scale; truncation scale factor that was used
688 */
689 {
690 int i, j;
691 int ny2, s10, s00, diff, dmax, dmin, s, smax;
692 int hm, h0, hp, hmm, hpm, hmp, hpp, hx2, hy2;
693 int m1,m2;
694
695 /*
696 * Maximum change in coefficients is determined by scale factor.
697 * Since we rounded during division (see digitize.c), the biggest
698 * permitted change is scale/2.
699 */
700 smax = (scale >> 1);
701 if (smax <= 0) return;
702 ny2 = ny << 1;
703 /*
704 * We're indexing a as a 2-D array with dimensions (nxtop,ny) of which
705 * only (nxtop,nytop) are used. The coefficients on the edge of the
706 * array are not adjusted (which is why the loops below start at 2
707 * instead of 0 and end at nxtop-2 instead of nxtop.)
708 */
709 /*
710 * Adjust x difference hx
711 */
712 for (i = 2; i<nxtop-2; i += 2) {
713 s00 = ny*i; /* s00 is index of a[i,j] */
714 s10 = s00+ny; /* s10 is index of a[i+1,j] */
715 for (j = 0; j<nytop; j += 2) {
716 /*
717 * hp is h0 (mean value) in next x zone, hm is h0 in previous x zone
718 */
719 hm = a[s00-ny2];
720 h0 = a[s00];
721 hp = a[s00+ny2];
722 /*
723 * diff = 8 * hx slope that would match h0 in neighboring zones
724 */
725 diff = hp-hm;
726 /*
727 * monotonicity constraints on diff
728 */
729 dmax = max( min( (hp-h0), (h0-hm) ), 0 ) << 2;
730 dmin = min( max( (hp-h0), (h0-hm) ), 0 ) << 2;
731 /*
732 * if monotonicity would set slope = 0 then don't change hx.
733 * note dmax>=0, dmin<=0.
734 */
735 if (dmin < dmax) {
736 diff = max( min(diff, dmax), dmin);
737 /*
738 * Compute change in slope limited to range +/- smax.
739 * Careful with rounding negative numbers when using
740 * shift for divide by 8.
741 */
742 s = diff-(a[s10]<<3);
743 s = (s>=0) ? (s>>3) : ((s+7)>>3) ;
744 s = max( min(s, smax), -smax);
745 a[s10] = a[s10]+s;
746 }
747 s00 += 2;
748 s10 += 2;
749 }
750 }
751 /*
752 * Adjust y difference hy
753 */
754 for (i = 0; i<nxtop; i += 2) {
755 s00 = ny*i+2;
756 s10 = s00+ny;
757 for (j = 2; j<nytop-2; j += 2) {
758 hm = a[s00-2];
759 h0 = a[s00];
760 hp = a[s00+2];
761 diff = hp-hm;
762 dmax = max( min( (hp-h0), (h0-hm) ), 0 ) << 2;
763 dmin = min( max( (hp-h0), (h0-hm) ), 0 ) << 2;
764 if (dmin < dmax) {
765 diff = max( min(diff, dmax), dmin);
766 s = diff-(a[s00+1]<<3);
767 s = (s>=0) ? (s>>3) : ((s+7)>>3) ;
768 s = max( min(s, smax), -smax);
769 a[s00+1] = a[s00+1]+s;
770 }
771 s00 += 2;
772 s10 += 2;
773 }
774 }
775 /*
776 * Adjust curvature difference hc
777 */
778 for (i = 2; i<nxtop-2; i += 2) {
779 s00 = ny*i+2;
780 s10 = s00+ny;
781 for (j = 2; j<nytop-2; j += 2) {
782 /*
783 * ------------------ y
784 * | hmp | | hpp | |
785 * ------------------ |
786 * | | h0 | | |
787 * ------------------ -------x
788 * | hmm | | hpm |
789 * ------------------
790 */
791 hmm = a[s00-ny2-2];
792 hpm = a[s00+ny2-2];
793 hmp = a[s00-ny2+2];
794 hpp = a[s00+ny2+2];
795 h0 = a[s00];
796 /*
797 * diff = 64 * hc value that would match h0 in neighboring zones
798 */
799 diff = hpp + hmm - hmp - hpm;
800 /*
801 * 2 times x,y slopes in this zone
802 */
803 hx2 = a[s10 ]<<1;
804 hy2 = a[s00+1]<<1;
805 /*
806 * monotonicity constraints on diff
807 */
808 m1 = min(max(hpp-h0,0)-hx2-hy2, max(h0-hpm,0)+hx2-hy2);
809 m2 = min(max(h0-hmp,0)-hx2+hy2, max(hmm-h0,0)+hx2+hy2);
810 dmax = min(m1,m2) << 4;
811 m1 = max(min(hpp-h0,0)-hx2-hy2, min(h0-hpm,0)+hx2-hy2);
812 m2 = max(min(h0-hmp,0)-hx2+hy2, min(hmm-h0,0)+hx2+hy2);
813 dmin = max(m1,m2) << 4;
814 /*
815 * if monotonicity would set slope = 0 then don't change hc.
816 * note dmax>=0, dmin<=0.
817 */
818 if (dmin < dmax) {
819 diff = max( min(diff, dmax), dmin);
820 /*
821 * Compute change in slope limited to range +/- smax.
822 * Careful with rounding negative numbers when using
823 * shift for divide by 64.
824 */
825 s = diff-(a[s10+1]<<6);
826 s = (s>=0) ? (s>>6) : ((s+63)>>6) ;
827 s = max( min(s, smax), -smax);
828 a[s10+1] = a[s10+1]+s;
829 }
830 s00 += 2;
831 s10 += 2;
832 }
833 }
834 }
835 /* ############################################################################ */
836 static void
hsmooth64(LONGLONG a[],int nxtop,int nytop,int ny,int scale)837 hsmooth64(LONGLONG a[], int nxtop, int nytop, int ny, int scale)
838 /*
839 LONGLONG a[]; array of H-transform coefficients
840 int nxtop,nytop; size of coefficient block to use
841 int ny; actual 1st dimension of array
842 int scale; truncation scale factor that was used
843 */
844 {
845 int i, j;
846 int ny2, s10, s00;
847 LONGLONG hm, h0, hp, hmm, hpm, hmp, hpp, hx2, hy2, diff, dmax, dmin, s, smax, m1, m2;
848
849 /*
850 * Maximum change in coefficients is determined by scale factor.
851 * Since we rounded during division (see digitize.c), the biggest
852 * permitted change is scale/2.
853 */
854 smax = (scale >> 1);
855 if (smax <= 0) return;
856 ny2 = ny << 1;
857 /*
858 * We're indexing a as a 2-D array with dimensions (nxtop,ny) of which
859 * only (nxtop,nytop) are used. The coefficients on the edge of the
860 * array are not adjusted (which is why the loops below start at 2
861 * instead of 0 and end at nxtop-2 instead of nxtop.)
862 */
863 /*
864 * Adjust x difference hx
865 */
866 for (i = 2; i<nxtop-2; i += 2) {
867 s00 = ny*i; /* s00 is index of a[i,j] */
868 s10 = s00+ny; /* s10 is index of a[i+1,j] */
869 for (j = 0; j<nytop; j += 2) {
870 /*
871 * hp is h0 (mean value) in next x zone, hm is h0 in previous x zone
872 */
873 hm = a[s00-ny2];
874 h0 = a[s00];
875 hp = a[s00+ny2];
876 /*
877 * diff = 8 * hx slope that would match h0 in neighboring zones
878 */
879 diff = hp-hm;
880 /*
881 * monotonicity constraints on diff
882 */
883 dmax = max( min( (hp-h0), (h0-hm) ), 0 ) << 2;
884 dmin = min( max( (hp-h0), (h0-hm) ), 0 ) << 2;
885 /*
886 * if monotonicity would set slope = 0 then don't change hx.
887 * note dmax>=0, dmin<=0.
888 */
889 if (dmin < dmax) {
890 diff = max( min(diff, dmax), dmin);
891 /*
892 * Compute change in slope limited to range +/- smax.
893 * Careful with rounding negative numbers when using
894 * shift for divide by 8.
895 */
896 s = diff-(a[s10]<<3);
897 s = (s>=0) ? (s>>3) : ((s+7)>>3) ;
898 s = max( min(s, smax), -smax);
899 a[s10] = a[s10]+s;
900 }
901 s00 += 2;
902 s10 += 2;
903 }
904 }
905 /*
906 * Adjust y difference hy
907 */
908 for (i = 0; i<nxtop; i += 2) {
909 s00 = ny*i+2;
910 s10 = s00+ny;
911 for (j = 2; j<nytop-2; j += 2) {
912 hm = a[s00-2];
913 h0 = a[s00];
914 hp = a[s00+2];
915 diff = hp-hm;
916 dmax = max( min( (hp-h0), (h0-hm) ), 0 ) << 2;
917 dmin = min( max( (hp-h0), (h0-hm) ), 0 ) << 2;
918 if (dmin < dmax) {
919 diff = max( min(diff, dmax), dmin);
920 s = diff-(a[s00+1]<<3);
921 s = (s>=0) ? (s>>3) : ((s+7)>>3) ;
922 s = max( min(s, smax), -smax);
923 a[s00+1] = a[s00+1]+s;
924 }
925 s00 += 2;
926 s10 += 2;
927 }
928 }
929 /*
930 * Adjust curvature difference hc
931 */
932 for (i = 2; i<nxtop-2; i += 2) {
933 s00 = ny*i+2;
934 s10 = s00+ny;
935 for (j = 2; j<nytop-2; j += 2) {
936 /*
937 * ------------------ y
938 * | hmp | | hpp | |
939 * ------------------ |
940 * | | h0 | | |
941 * ------------------ -------x
942 * | hmm | | hpm |
943 * ------------------
944 */
945 hmm = a[s00-ny2-2];
946 hpm = a[s00+ny2-2];
947 hmp = a[s00-ny2+2];
948 hpp = a[s00+ny2+2];
949 h0 = a[s00];
950 /*
951 * diff = 64 * hc value that would match h0 in neighboring zones
952 */
953 diff = hpp + hmm - hmp - hpm;
954 /*
955 * 2 times x,y slopes in this zone
956 */
957 hx2 = a[s10 ]<<1;
958 hy2 = a[s00+1]<<1;
959 /*
960 * monotonicity constraints on diff
961 */
962 m1 = min(max(hpp-h0,0)-hx2-hy2, max(h0-hpm,0)+hx2-hy2);
963 m2 = min(max(h0-hmp,0)-hx2+hy2, max(hmm-h0,0)+hx2+hy2);
964 dmax = min(m1,m2) << 4;
965 m1 = max(min(hpp-h0,0)-hx2-hy2, min(h0-hpm,0)+hx2-hy2);
966 m2 = max(min(h0-hmp,0)-hx2+hy2, min(hmm-h0,0)+hx2+hy2);
967 dmin = max(m1,m2) << 4;
968 /*
969 * if monotonicity would set slope = 0 then don't change hc.
970 * note dmax>=0, dmin<=0.
971 */
972 if (dmin < dmax) {
973 diff = max( min(diff, dmax), dmin);
974 /*
975 * Compute change in slope limited to range +/- smax.
976 * Careful with rounding negative numbers when using
977 * shift for divide by 64.
978 */
979 s = diff-(a[s10+1]<<6);
980 s = (s>=0) ? (s>>6) : ((s+63)>>6) ;
981 s = max( min(s, smax), -smax);
982 a[s10+1] = a[s10+1]+s;
983 }
984 s00 += 2;
985 s10 += 2;
986 }
987 }
988 }
989
990
991 /* ############################################################################ */
992 /* ############################################################################ */
993 /* Copyright (c) 1993 Association of Universities for Research
994 * in Astronomy. All rights reserved. Produced under National
995 * Aeronautics and Space Administration Contract No. NAS5-26555.
996 */
997 /* undigitize.c undigitize H-transform
998 *
999 * Programmer: R. White Date: 9 May 1991
1000 */
1001
1002 /* ############################################################################ */
1003 static void
undigitize(int a[],int nx,int ny,int scale)1004 undigitize(int a[], int nx, int ny, int scale)
1005 {
1006 int *p;
1007
1008 /*
1009 * multiply by scale
1010 */
1011 if (scale <= 1) return;
1012 for (p=a; p <= &a[nx*ny-1]; p++) *p = (*p)*scale;
1013 }
1014 /* ############################################################################ */
1015 static void
undigitize64(LONGLONG a[],int nx,int ny,int scale)1016 undigitize64(LONGLONG a[], int nx, int ny, int scale)
1017 {
1018 LONGLONG *p, scale64;
1019
1020 /*
1021 * multiply by scale
1022 */
1023 if (scale <= 1) return;
1024 scale64 = (LONGLONG) scale; /* use a 64-bit int for efficiency in the big loop */
1025
1026 for (p=a; p <= &a[nx*ny-1]; p++) *p = (*p)*scale64;
1027 }
1028
1029 /* ############################################################################ */
1030 /* ############################################################################ */
1031 /* Copyright (c) 1993 Association of Universities for Research
1032 * in Astronomy. All rights reserved. Produced under National
1033 * Aeronautics and Space Administration Contract No. NAS5-26555.
1034 */
1035 /* decode.c read codes from infile and construct array
1036 *
1037 * Programmer: R. White Date: 2 February 1994
1038 */
1039
1040
1041 static char code_magic[2] = { (char)0xDD, (char)0x99 };
1042
1043 /* ############################################################################ */
decode(unsigned char * infile,int * a,int * nx,int * ny,int * scale)1044 static int decode(unsigned char *infile, int *a, int *nx, int *ny, int *scale)
1045 /*
1046 char *infile; input file
1047 int *a; address of output array [nx][ny]
1048 int *nx,*ny; size of output array
1049 int *scale; scale factor for digitization
1050 */
1051 {
1052 LONGLONG sumall;
1053 int stat;
1054 unsigned char nbitplanes[3];
1055 char tmagic[2];
1056
1057 /* initialize the byte read position to the beginning of the array */;
1058 nextchar = 0;
1059
1060 /*
1061 * File starts either with special 2-byte magic code or with
1062 * FITS keyword "SIMPLE ="
1063 */
1064 qread(infile, tmagic, sizeof(tmagic));
1065 /*
1066 * check for correct magic code value
1067 */
1068 if (memcmp(tmagic,code_magic,sizeof(code_magic)) != 0) {
1069 ffpmsg("bad file format");
1070 return(DATA_DECOMPRESSION_ERR);
1071 }
1072 *nx =readint(infile); /* x size of image */
1073 *ny =readint(infile); /* y size of image */
1074 *scale=readint(infile); /* scale factor for digitization */
1075
1076 /* sum of all pixels */
1077 sumall=readlonglong(infile);
1078 /* # bits in quadrants */
1079
1080 qread(infile, (char *) nbitplanes, sizeof(nbitplanes));
1081
1082 stat = dodecode(infile, a, *nx, *ny, nbitplanes);
1083 /*
1084 * put sum of all pixels back into pixel 0
1085 */
1086 a[0] = (int) sumall;
1087 return(stat);
1088 }
1089 /* ############################################################################ */
decode64(unsigned char * infile,LONGLONG * a,int * nx,int * ny,int * scale)1090 static int decode64(unsigned char *infile, LONGLONG *a, int *nx, int *ny, int *scale)
1091 /*
1092 char *infile; input file
1093 LONGLONG *a; address of output array [nx][ny]
1094 int *nx,*ny; size of output array
1095 int *scale; scale factor for digitization
1096 */
1097 {
1098 int stat;
1099 LONGLONG sumall;
1100 unsigned char nbitplanes[3];
1101 char tmagic[2];
1102
1103 /* initialize the byte read position to the beginning of the array */;
1104 nextchar = 0;
1105
1106 /*
1107 * File starts either with special 2-byte magic code or with
1108 * FITS keyword "SIMPLE ="
1109 */
1110 qread(infile, tmagic, sizeof(tmagic));
1111 /*
1112 * check for correct magic code value
1113 */
1114 if (memcmp(tmagic,code_magic,sizeof(code_magic)) != 0) {
1115 ffpmsg("bad file format");
1116 return(DATA_DECOMPRESSION_ERR);
1117 }
1118 *nx =readint(infile); /* x size of image */
1119 *ny =readint(infile); /* y size of image */
1120 *scale=readint(infile); /* scale factor for digitization */
1121
1122 /* sum of all pixels */
1123 sumall=readlonglong(infile);
1124 /* # bits in quadrants */
1125
1126 qread(infile, (char *) nbitplanes, sizeof(nbitplanes));
1127
1128 stat = dodecode64(infile, a, *nx, *ny, nbitplanes);
1129 /*
1130 * put sum of all pixels back into pixel 0
1131 */
1132 a[0] = sumall;
1133
1134 return(stat);
1135 }
1136
1137
1138 /* ############################################################################ */
1139 /* ############################################################################ */
1140 /* Copyright (c) 1993 Association of Universities for Research
1141 * in Astronomy. All rights reserved. Produced under National
1142 * Aeronautics and Space Administration Contract No. NAS5-26555.
1143 */
1144 /* dodecode.c Decode stream of characters on infile and return array
1145 *
1146 * This version encodes the different quadrants separately
1147 *
1148 * Programmer: R. White Date: 9 May 1991
1149 */
1150
1151 /* ############################################################################ */
1152 static int
dodecode(unsigned char * infile,int a[],int nx,int ny,unsigned char nbitplanes[3])1153 dodecode(unsigned char *infile, int a[], int nx, int ny, unsigned char nbitplanes[3])
1154
1155 /* int a[];
1156 int nx,ny; Array dimensions are [nx][ny]
1157 unsigned char nbitplanes[3]; Number of bit planes in quadrants
1158 */
1159 {
1160 int i, nel, nx2, ny2, stat;
1161
1162 nel = nx*ny;
1163 nx2 = (nx+1)/2;
1164 ny2 = (ny+1)/2;
1165
1166 /*
1167 * initialize a to zero
1168 */
1169 for (i=0; i<nel; i++) a[i] = 0;
1170 /*
1171 * Initialize bit input
1172 */
1173 start_inputing_bits();
1174 /*
1175 * read bit planes for each quadrant
1176 */
1177 stat = qtree_decode(infile, &a[0], ny, nx2, ny2, nbitplanes[0]);
1178 if (stat) return(stat);
1179
1180 stat = qtree_decode(infile, &a[ny2], ny, nx2, ny/2, nbitplanes[1]);
1181 if (stat) return(stat);
1182
1183 stat = qtree_decode(infile, &a[ny*nx2], ny, nx/2, ny2, nbitplanes[1]);
1184 if (stat) return(stat);
1185
1186 stat = qtree_decode(infile, &a[ny*nx2+ny2], ny, nx/2, ny/2, nbitplanes[2]);
1187 if (stat) return(stat);
1188
1189 /*
1190 * make sure there is an EOF symbol (nybble=0) at end
1191 */
1192 if (input_nybble(infile) != 0) {
1193 ffpmsg("dodecode: bad bit plane values");
1194 return(DATA_DECOMPRESSION_ERR);
1195 }
1196 /*
1197 * now get the sign bits
1198 * Re-initialize bit input
1199 */
1200 start_inputing_bits();
1201 for (i=0; i<nel; i++) {
1202 if (a[i]) {
1203 /* tried putting the input_bit code in-line here, instead of */
1204 /* calling the function, but it made no difference in the speed */
1205 if (input_bit(infile)) a[i] = -a[i];
1206 }
1207 }
1208 return(0);
1209 }
1210 /* ############################################################################ */
1211 static int
dodecode64(unsigned char * infile,LONGLONG a[],int nx,int ny,unsigned char nbitplanes[3])1212 dodecode64(unsigned char *infile, LONGLONG a[], int nx, int ny, unsigned char nbitplanes[3])
1213
1214 /* LONGLONG a[];
1215 int nx,ny; Array dimensions are [nx][ny]
1216 unsigned char nbitplanes[3]; Number of bit planes in quadrants
1217 */
1218 {
1219 int i, nel, nx2, ny2, stat;
1220
1221 nel = nx*ny;
1222 nx2 = (nx+1)/2;
1223 ny2 = (ny+1)/2;
1224
1225 /*
1226 * initialize a to zero
1227 */
1228 for (i=0; i<nel; i++) a[i] = 0;
1229 /*
1230 * Initialize bit input
1231 */
1232 start_inputing_bits();
1233 /*
1234 * read bit planes for each quadrant
1235 */
1236 stat = qtree_decode64(infile, &a[0], ny, nx2, ny2, nbitplanes[0]);
1237 if (stat) return(stat);
1238
1239 stat = qtree_decode64(infile, &a[ny2], ny, nx2, ny/2, nbitplanes[1]);
1240 if (stat) return(stat);
1241
1242 stat = qtree_decode64(infile, &a[ny*nx2], ny, nx/2, ny2, nbitplanes[1]);
1243 if (stat) return(stat);
1244
1245 stat = qtree_decode64(infile, &a[ny*nx2+ny2], ny, nx/2, ny/2, nbitplanes[2]);
1246 if (stat) return(stat);
1247
1248 /*
1249 * make sure there is an EOF symbol (nybble=0) at end
1250 */
1251 if (input_nybble(infile) != 0) {
1252 ffpmsg("dodecode64: bad bit plane values");
1253 return(DATA_DECOMPRESSION_ERR);
1254 }
1255 /*
1256 * now get the sign bits
1257 * Re-initialize bit input
1258 */
1259 start_inputing_bits();
1260 for (i=0; i<nel; i++) {
1261 if (a[i]) {
1262 if (input_bit(infile) != 0) a[i] = -a[i];
1263 }
1264 }
1265 return(0);
1266 }
1267
1268 /* ############################################################################ */
1269 /* ############################################################################ */
1270 /* Copyright (c) 1993 Association of Universities for Research
1271 * in Astronomy. All rights reserved. Produced under National
1272 * Aeronautics and Space Administration Contract No. NAS5-26555.
1273 */
1274 /* qtree_decode.c Read stream of codes from infile and construct bit planes
1275 * in quadrant of 2-D array using binary quadtree coding
1276 *
1277 * Programmer: R. White Date: 7 May 1991
1278 */
1279
1280 /* ############################################################################ */
1281 static int
qtree_decode(unsigned char * infile,int a[],int n,int nqx,int nqy,int nbitplanes)1282 qtree_decode(unsigned char *infile, int a[], int n, int nqx, int nqy, int nbitplanes)
1283
1284 /*
1285 char *infile;
1286 int a[]; a is 2-D array with dimensions (n,n)
1287 int n; length of full row in a
1288 int nqx; partial length of row to decode
1289 int nqy; partial length of column (<=n)
1290 int nbitplanes; number of bitplanes to decode
1291 */
1292 {
1293 int log2n, k, bit, b, nqmax;
1294 int nx,ny,nfx,nfy,c;
1295 int nqx2, nqy2;
1296 unsigned char *scratch;
1297
1298 /*
1299 * log2n is log2 of max(nqx,nqy) rounded up to next power of 2
1300 */
1301 nqmax = (nqx>nqy) ? nqx : nqy;
1302 log2n = (int) (log((float) nqmax)/log(2.0)+0.5);
1303 if (nqmax > (1<<log2n)) {
1304 log2n += 1;
1305 }
1306 /*
1307 * allocate scratch array for working space
1308 */
1309 nqx2=(nqx+1)/2;
1310 nqy2=(nqy+1)/2;
1311 scratch = (unsigned char *) malloc(nqx2*nqy2);
1312 if (scratch == (unsigned char *) NULL) {
1313 ffpmsg("qtree_decode: insufficient memory");
1314 return(DATA_DECOMPRESSION_ERR);
1315 }
1316 /*
1317 * now decode each bit plane, starting at the top
1318 * A is assumed to be initialized to zero
1319 */
1320 for (bit = nbitplanes-1; bit >= 0; bit--) {
1321 /*
1322 * Was bitplane was quadtree-coded or written directly?
1323 */
1324 b = input_nybble(infile);
1325
1326 if(b == 0) {
1327 /*
1328 * bit map was written directly
1329 */
1330 read_bdirect(infile,a,n,nqx,nqy,scratch,bit);
1331 } else if (b != 0xf) {
1332 ffpmsg("qtree_decode: bad format code");
1333 return(DATA_DECOMPRESSION_ERR);
1334 } else {
1335 /*
1336 * bitmap was quadtree-coded, do log2n expansions
1337 *
1338 * read first code
1339 */
1340 scratch[0] = input_huffman(infile);
1341 /*
1342 * now do log2n expansions, reading codes from file as necessary
1343 */
1344 nx = 1;
1345 ny = 1;
1346 nfx = nqx;
1347 nfy = nqy;
1348 c = 1<<log2n;
1349 for (k = 1; k<log2n; k++) {
1350 /*
1351 * this somewhat cryptic code generates the sequence
1352 * n[k-1] = (n[k]+1)/2 where n[log2n]=nqx or nqy
1353 */
1354 c = c>>1;
1355 nx = nx<<1;
1356 ny = ny<<1;
1357 if (nfx <= c) { nx -= 1; } else { nfx -= c; }
1358 if (nfy <= c) { ny -= 1; } else { nfy -= c; }
1359 qtree_expand(infile,scratch,nx,ny,scratch);
1360 }
1361 /*
1362 * now copy last set of 4-bit codes to bitplane bit of array a
1363 */
1364 qtree_bitins(scratch,nqx,nqy,a,n,bit);
1365 }
1366 }
1367 free(scratch);
1368 return(0);
1369 }
1370 /* ############################################################################ */
1371 static int
qtree_decode64(unsigned char * infile,LONGLONG a[],int n,int nqx,int nqy,int nbitplanes)1372 qtree_decode64(unsigned char *infile, LONGLONG a[], int n, int nqx, int nqy, int nbitplanes)
1373
1374 /*
1375 char *infile;
1376 LONGLONG a[]; a is 2-D array with dimensions (n,n)
1377 int n; length of full row in a
1378 int nqx; partial length of row to decode
1379 int nqy; partial length of column (<=n)
1380 int nbitplanes; number of bitplanes to decode
1381 */
1382 {
1383 int log2n, k, bit, b, nqmax;
1384 int nx,ny,nfx,nfy,c;
1385 int nqx2, nqy2;
1386 unsigned char *scratch;
1387
1388 /*
1389 * log2n is log2 of max(nqx,nqy) rounded up to next power of 2
1390 */
1391 nqmax = (nqx>nqy) ? nqx : nqy;
1392 log2n = (int) (log((float) nqmax)/log(2.0)+0.5);
1393 if (nqmax > (1<<log2n)) {
1394 log2n += 1;
1395 }
1396 /*
1397 * allocate scratch array for working space
1398 */
1399 nqx2=(nqx+1)/2;
1400 nqy2=(nqy+1)/2;
1401 scratch = (unsigned char *) malloc(nqx2*nqy2);
1402 if (scratch == (unsigned char *) NULL) {
1403 ffpmsg("qtree_decode64: insufficient memory");
1404 return(DATA_DECOMPRESSION_ERR);
1405 }
1406 /*
1407 * now decode each bit plane, starting at the top
1408 * A is assumed to be initialized to zero
1409 */
1410 for (bit = nbitplanes-1; bit >= 0; bit--) {
1411 /*
1412 * Was bitplane was quadtree-coded or written directly?
1413 */
1414 b = input_nybble(infile);
1415
1416 if(b == 0) {
1417 /*
1418 * bit map was written directly
1419 */
1420 read_bdirect64(infile,a,n,nqx,nqy,scratch,bit);
1421 } else if (b != 0xf) {
1422 ffpmsg("qtree_decode64: bad format code");
1423 return(DATA_DECOMPRESSION_ERR);
1424 } else {
1425 /*
1426 * bitmap was quadtree-coded, do log2n expansions
1427 *
1428 * read first code
1429 */
1430 scratch[0] = input_huffman(infile);
1431 /*
1432 * now do log2n expansions, reading codes from file as necessary
1433 */
1434 nx = 1;
1435 ny = 1;
1436 nfx = nqx;
1437 nfy = nqy;
1438 c = 1<<log2n;
1439 for (k = 1; k<log2n; k++) {
1440 /*
1441 * this somewhat cryptic code generates the sequence
1442 * n[k-1] = (n[k]+1)/2 where n[log2n]=nqx or nqy
1443 */
1444 c = c>>1;
1445 nx = nx<<1;
1446 ny = ny<<1;
1447 if (nfx <= c) { nx -= 1; } else { nfx -= c; }
1448 if (nfy <= c) { ny -= 1; } else { nfy -= c; }
1449 qtree_expand(infile,scratch,nx,ny,scratch);
1450 }
1451 /*
1452 * now copy last set of 4-bit codes to bitplane bit of array a
1453 */
1454 qtree_bitins64(scratch,nqx,nqy,a,n,bit);
1455 }
1456 }
1457 free(scratch);
1458 return(0);
1459 }
1460
1461
1462 /* ############################################################################ */
1463 /*
1464 * do one quadtree expansion step on array a[(nqx+1)/2,(nqy+1)/2]
1465 * results put into b[nqx,nqy] (which may be the same as a)
1466 */
1467 static void
qtree_expand(unsigned char * infile,unsigned char a[],int nx,int ny,unsigned char b[])1468 qtree_expand(unsigned char *infile, unsigned char a[], int nx, int ny, unsigned char b[])
1469 {
1470 int i;
1471
1472 /*
1473 * first copy a to b, expanding each 4-bit value
1474 */
1475 qtree_copy(a,nx,ny,b,ny);
1476 /*
1477 * now read new 4-bit values into b for each non-zero element
1478 */
1479 for (i = nx*ny-1; i >= 0; i--) {
1480 if (b[i]) b[i] = input_huffman(infile);
1481 }
1482 }
1483
1484 /* ############################################################################ */
1485 /*
1486 * copy 4-bit values from a[(nx+1)/2,(ny+1)/2] to b[nx,ny], expanding
1487 * each value to 2x2 pixels
1488 * a,b may be same array
1489 */
1490 static void
qtree_copy(unsigned char a[],int nx,int ny,unsigned char b[],int n)1491 qtree_copy(unsigned char a[], int nx, int ny, unsigned char b[], int n)
1492 /* int n; declared y dimension of b */
1493 {
1494 int i, j, k, nx2, ny2;
1495 int s00, s10;
1496
1497 /*
1498 * first copy 4-bit values to b
1499 * start at end in case a,b are same array
1500 */
1501 nx2 = (nx+1)/2;
1502 ny2 = (ny+1)/2;
1503 k = ny2*(nx2-1)+ny2-1; /* k is index of a[i,j] */
1504 for (i = nx2-1; i >= 0; i--) {
1505 s00 = 2*(n*i+ny2-1); /* s00 is index of b[2*i,2*j] */
1506 for (j = ny2-1; j >= 0; j--) {
1507 b[s00] = a[k];
1508 k -= 1;
1509 s00 -= 2;
1510 }
1511 }
1512 /*
1513 * now expand each 2x2 block
1514 */
1515 for (i = 0; i<nx-1; i += 2) {
1516
1517 /* Note:
1518 Unlike the case in qtree_bitins, this code runs faster on a 32-bit linux
1519 machine using the s10 intermediate variable, rather that using s00+n.
1520 Go figure!
1521 */
1522 s00 = n*i; /* s00 is index of b[i,j] */
1523 s10 = s00+n; /* s10 is index of b[i+1,j] */
1524
1525 for (j = 0; j<ny-1; j += 2) {
1526
1527 switch (b[s00]) {
1528 case(0):
1529 b[s10+1] = 0;
1530 b[s10 ] = 0;
1531 b[s00+1] = 0;
1532 b[s00 ] = 0;
1533
1534 break;
1535 case(1):
1536 b[s10+1] = 1;
1537 b[s10 ] = 0;
1538 b[s00+1] = 0;
1539 b[s00 ] = 0;
1540
1541 break;
1542 case(2):
1543 b[s10+1] = 0;
1544 b[s10 ] = 1;
1545 b[s00+1] = 0;
1546 b[s00 ] = 0;
1547
1548 break;
1549 case(3):
1550 b[s10+1] = 1;
1551 b[s10 ] = 1;
1552 b[s00+1] = 0;
1553 b[s00 ] = 0;
1554
1555 break;
1556 case(4):
1557 b[s10+1] = 0;
1558 b[s10 ] = 0;
1559 b[s00+1] = 1;
1560 b[s00 ] = 0;
1561
1562 break;
1563 case(5):
1564 b[s10+1] = 1;
1565 b[s10 ] = 0;
1566 b[s00+1] = 1;
1567 b[s00 ] = 0;
1568
1569 break;
1570 case(6):
1571 b[s10+1] = 0;
1572 b[s10 ] = 1;
1573 b[s00+1] = 1;
1574 b[s00 ] = 0;
1575
1576 break;
1577 case(7):
1578 b[s10+1] = 1;
1579 b[s10 ] = 1;
1580 b[s00+1] = 1;
1581 b[s00 ] = 0;
1582
1583 break;
1584 case(8):
1585 b[s10+1] = 0;
1586 b[s10 ] = 0;
1587 b[s00+1] = 0;
1588 b[s00 ] = 1;
1589
1590 break;
1591 case(9):
1592 b[s10+1] = 1;
1593 b[s10 ] = 0;
1594 b[s00+1] = 0;
1595 b[s00 ] = 1;
1596 break;
1597 case(10):
1598 b[s10+1] = 0;
1599 b[s10 ] = 1;
1600 b[s00+1] = 0;
1601 b[s00 ] = 1;
1602
1603 break;
1604 case(11):
1605 b[s10+1] = 1;
1606 b[s10 ] = 1;
1607 b[s00+1] = 0;
1608 b[s00 ] = 1;
1609
1610 break;
1611 case(12):
1612 b[s10+1] = 0;
1613 b[s10 ] = 0;
1614 b[s00+1] = 1;
1615 b[s00 ] = 1;
1616
1617 break;
1618 case(13):
1619 b[s10+1] = 1;
1620 b[s10 ] = 0;
1621 b[s00+1] = 1;
1622 b[s00 ] = 1;
1623
1624 break;
1625 case(14):
1626 b[s10+1] = 0;
1627 b[s10 ] = 1;
1628 b[s00+1] = 1;
1629 b[s00 ] = 1;
1630
1631 break;
1632 case(15):
1633 b[s10+1] = 1;
1634 b[s10 ] = 1;
1635 b[s00+1] = 1;
1636 b[s00 ] = 1;
1637
1638 break;
1639 }
1640 /*
1641 b[s10+1] = b[s00] & 1;
1642 b[s10 ] = (b[s00]>>1) & 1;
1643 b[s00+1] = (b[s00]>>2) & 1;
1644 b[s00 ] = (b[s00]>>3) & 1;
1645 */
1646
1647 s00 += 2;
1648 s10 += 2;
1649 }
1650
1651 if (j < ny) {
1652 /*
1653 * row size is odd, do last element in row
1654 * s00+1, s10+1 are off edge
1655 */
1656 /* not worth converting this to use 16 case statements */
1657 b[s10 ] = (b[s00]>>1) & 1;
1658 b[s00 ] = (b[s00]>>3) & 1;
1659 }
1660 }
1661 if (i < nx) {
1662 /*
1663 * column size is odd, do last row
1664 * s10, s10+1 are off edge
1665 */
1666 s00 = n*i;
1667 for (j = 0; j<ny-1; j += 2) {
1668 /* not worth converting this to use 16 case statements */
1669 b[s00+1] = (b[s00]>>2) & 1;
1670 b[s00 ] = (b[s00]>>3) & 1;
1671 s00 += 2;
1672 }
1673 if (j < ny) {
1674 /*
1675 * both row and column size are odd, do corner element
1676 * s00+1, s10, s10+1 are off edge
1677 */
1678 /* not worth converting this to use 16 case statements */
1679 b[s00 ] = (b[s00]>>3) & 1;
1680 }
1681 }
1682 }
1683
1684 /* ############################################################################ */
1685 /*
1686 * Copy 4-bit values from a[(nx+1)/2,(ny+1)/2] to b[nx,ny], expanding
1687 * each value to 2x2 pixels and inserting into bitplane BIT of B.
1688 * A,B may NOT be same array (it wouldn't make sense to be inserting
1689 * bits into the same array anyway.)
1690 */
1691 static void
qtree_bitins(unsigned char a[],int nx,int ny,int b[],int n,int bit)1692 qtree_bitins(unsigned char a[], int nx, int ny, int b[], int n, int bit)
1693 /*
1694 int n; declared y dimension of b
1695 */
1696 {
1697 int i, j, k;
1698 int s00;
1699 int plane_val;
1700
1701 plane_val = 1 << bit;
1702
1703 /*
1704 * expand each 2x2 block
1705 */
1706 k = 0; /* k is index of a[i/2,j/2] */
1707 for (i = 0; i<nx-1; i += 2) {
1708 s00 = n*i; /* s00 is index of b[i,j] */
1709
1710 /* Note:
1711 this code appears to run very slightly faster on a 32-bit linux
1712 machine using s00+n rather than the s10 intermediate variable
1713 */
1714 /* s10 = s00+n; */ /* s10 is index of b[i+1,j] */
1715 for (j = 0; j<ny-1; j += 2) {
1716
1717 switch (a[k]) {
1718 case(0):
1719 break;
1720 case(1):
1721 b[s00+n+1] |= plane_val;
1722 break;
1723 case(2):
1724 b[s00+n ] |= plane_val;
1725 break;
1726 case(3):
1727 b[s00+n+1] |= plane_val;
1728 b[s00+n ] |= plane_val;
1729 break;
1730 case(4):
1731 b[s00+1] |= plane_val;
1732 break;
1733 case(5):
1734 b[s00+n+1] |= plane_val;
1735 b[s00+1] |= plane_val;
1736 break;
1737 case(6):
1738 b[s00+n ] |= plane_val;
1739 b[s00+1] |= plane_val;
1740 break;
1741 case(7):
1742 b[s00+n+1] |= plane_val;
1743 b[s00+n ] |= plane_val;
1744 b[s00+1] |= plane_val;
1745 break;
1746 case(8):
1747 b[s00 ] |= plane_val;
1748 break;
1749 case(9):
1750 b[s00+n+1] |= plane_val;
1751 b[s00 ] |= plane_val;
1752 break;
1753 case(10):
1754 b[s00+n ] |= plane_val;
1755 b[s00 ] |= plane_val;
1756 break;
1757 case(11):
1758 b[s00+n+1] |= plane_val;
1759 b[s00+n ] |= plane_val;
1760 b[s00 ] |= plane_val;
1761 break;
1762 case(12):
1763 b[s00+1] |= plane_val;
1764 b[s00 ] |= plane_val;
1765 break;
1766 case(13):
1767 b[s00+n+1] |= plane_val;
1768 b[s00+1] |= plane_val;
1769 b[s00 ] |= plane_val;
1770 break;
1771 case(14):
1772 b[s00+n ] |= plane_val;
1773 b[s00+1] |= plane_val;
1774 b[s00 ] |= plane_val;
1775 break;
1776 case(15):
1777 b[s00+n+1] |= plane_val;
1778 b[s00+n ] |= plane_val;
1779 b[s00+1] |= plane_val;
1780 b[s00 ] |= plane_val;
1781 break;
1782 }
1783
1784 /*
1785 b[s10+1] |= ( a[k] & 1) << bit;
1786 b[s10 ] |= ((a[k]>>1) & 1) << bit;
1787 b[s00+1] |= ((a[k]>>2) & 1) << bit;
1788 b[s00 ] |= ((a[k]>>3) & 1) << bit;
1789 */
1790 s00 += 2;
1791 /* s10 += 2; */
1792 k += 1;
1793 }
1794 if (j < ny) {
1795 /*
1796 * row size is odd, do last element in row
1797 * s00+1, s10+1 are off edge
1798 */
1799
1800 switch (a[k]) {
1801 case(0):
1802 break;
1803 case(1):
1804 break;
1805 case(2):
1806 b[s00+n ] |= plane_val;
1807 break;
1808 case(3):
1809 b[s00+n ] |= plane_val;
1810 break;
1811 case(4):
1812 break;
1813 case(5):
1814 break;
1815 case(6):
1816 b[s00+n ] |= plane_val;
1817 break;
1818 case(7):
1819 b[s00+n ] |= plane_val;
1820 break;
1821 case(8):
1822 b[s00 ] |= plane_val;
1823 break;
1824 case(9):
1825 b[s00 ] |= plane_val;
1826 break;
1827 case(10):
1828 b[s00+n ] |= plane_val;
1829 b[s00 ] |= plane_val;
1830 break;
1831 case(11):
1832 b[s00+n ] |= plane_val;
1833 b[s00 ] |= plane_val;
1834 break;
1835 case(12):
1836 b[s00 ] |= plane_val;
1837 break;
1838 case(13):
1839 b[s00 ] |= plane_val;
1840 break;
1841 case(14):
1842 b[s00+n ] |= plane_val;
1843 b[s00 ] |= plane_val;
1844 break;
1845 case(15):
1846 b[s00+n ] |= plane_val;
1847 b[s00 ] |= plane_val;
1848 break;
1849 }
1850
1851 /*
1852 b[s10 ] |= ((a[k]>>1) & 1) << bit;
1853 b[s00 ] |= ((a[k]>>3) & 1) << bit;
1854 */
1855 k += 1;
1856 }
1857 }
1858 if (i < nx) {
1859 /*
1860 * column size is odd, do last row
1861 * s10, s10+1 are off edge
1862 */
1863 s00 = n*i;
1864 for (j = 0; j<ny-1; j += 2) {
1865
1866 switch (a[k]) {
1867 case(0):
1868 break;
1869 case(1):
1870 break;
1871 case(2):
1872 break;
1873 case(3):
1874 break;
1875 case(4):
1876 b[s00+1] |= plane_val;
1877 break;
1878 case(5):
1879 b[s00+1] |= plane_val;
1880 break;
1881 case(6):
1882 b[s00+1] |= plane_val;
1883 break;
1884 case(7):
1885 b[s00+1] |= plane_val;
1886 break;
1887 case(8):
1888 b[s00 ] |= plane_val;
1889 break;
1890 case(9):
1891 b[s00 ] |= plane_val;
1892 break;
1893 case(10):
1894 b[s00 ] |= plane_val;
1895 break;
1896 case(11):
1897 b[s00 ] |= plane_val;
1898 break;
1899 case(12):
1900 b[s00+1] |= plane_val;
1901 b[s00 ] |= plane_val;
1902 break;
1903 case(13):
1904 b[s00+1] |= plane_val;
1905 b[s00 ] |= plane_val;
1906 break;
1907 case(14):
1908 b[s00+1] |= plane_val;
1909 b[s00 ] |= plane_val;
1910 break;
1911 case(15):
1912 b[s00+1] |= plane_val;
1913 b[s00 ] |= plane_val;
1914 break;
1915 }
1916
1917 /*
1918 b[s00+1] |= ((a[k]>>2) & 1) << bit;
1919 b[s00 ] |= ((a[k]>>3) & 1) << bit;
1920 */
1921
1922 s00 += 2;
1923 k += 1;
1924 }
1925 if (j < ny) {
1926 /*
1927 * both row and column size are odd, do corner element
1928 * s00+1, s10, s10+1 are off edge
1929 */
1930
1931 switch (a[k]) {
1932 case(0):
1933 break;
1934 case(1):
1935 break;
1936 case(2):
1937 break;
1938 case(3):
1939 break;
1940 case(4):
1941 break;
1942 case(5):
1943 break;
1944 case(6):
1945 break;
1946 case(7):
1947 break;
1948 case(8):
1949 b[s00 ] |= plane_val;
1950 break;
1951 case(9):
1952 b[s00 ] |= plane_val;
1953 break;
1954 case(10):
1955 b[s00 ] |= plane_val;
1956 break;
1957 case(11):
1958 b[s00 ] |= plane_val;
1959 break;
1960 case(12):
1961 b[s00 ] |= plane_val;
1962 break;
1963 case(13):
1964 b[s00 ] |= plane_val;
1965 break;
1966 case(14):
1967 b[s00 ] |= plane_val;
1968 break;
1969 case(15):
1970 b[s00 ] |= plane_val;
1971 break;
1972 }
1973
1974 /*
1975 b[s00 ] |= ((a[k]>>3) & 1) << bit;
1976 */
1977 k += 1;
1978 }
1979 }
1980 }
1981 /* ############################################################################ */
1982 /*
1983 * Copy 4-bit values from a[(nx+1)/2,(ny+1)/2] to b[nx,ny], expanding
1984 * each value to 2x2 pixels and inserting into bitplane BIT of B.
1985 * A,B may NOT be same array (it wouldn't make sense to be inserting
1986 * bits into the same array anyway.)
1987 */
1988 static void
qtree_bitins64(unsigned char a[],int nx,int ny,LONGLONG b[],int n,int bit)1989 qtree_bitins64(unsigned char a[], int nx, int ny, LONGLONG b[], int n, int bit)
1990 /*
1991 int n; declared y dimension of b
1992 */
1993 {
1994 int i, j, k;
1995 int s00;
1996 LONGLONG plane_val;
1997
1998 plane_val = ((LONGLONG) 1) << bit;
1999
2000 /*
2001 * expand each 2x2 block
2002 */
2003 k = 0; /* k is index of a[i/2,j/2] */
2004 for (i = 0; i<nx-1; i += 2) {
2005 s00 = n*i; /* s00 is index of b[i,j] */
2006
2007 /* Note:
2008 this code appears to run very slightly faster on a 32-bit linux
2009 machine using s00+n rather than the s10 intermediate variable
2010 */
2011 /* s10 = s00+n; */ /* s10 is index of b[i+1,j] */
2012 for (j = 0; j<ny-1; j += 2) {
2013
2014 switch (a[k]) {
2015 case(0):
2016 break;
2017 case(1):
2018 b[s00+n+1] |= plane_val;
2019 break;
2020 case(2):
2021 b[s00+n ] |= plane_val;
2022 break;
2023 case(3):
2024 b[s00+n+1] |= plane_val;
2025 b[s00+n ] |= plane_val;
2026 break;
2027 case(4):
2028 b[s00+1] |= plane_val;
2029 break;
2030 case(5):
2031 b[s00+n+1] |= plane_val;
2032 b[s00+1] |= plane_val;
2033 break;
2034 case(6):
2035 b[s00+n ] |= plane_val;
2036 b[s00+1] |= plane_val;
2037 break;
2038 case(7):
2039 b[s00+n+1] |= plane_val;
2040 b[s00+n ] |= plane_val;
2041 b[s00+1] |= plane_val;
2042 break;
2043 case(8):
2044 b[s00 ] |= plane_val;
2045 break;
2046 case(9):
2047 b[s00+n+1] |= plane_val;
2048 b[s00 ] |= plane_val;
2049 break;
2050 case(10):
2051 b[s00+n ] |= plane_val;
2052 b[s00 ] |= plane_val;
2053 break;
2054 case(11):
2055 b[s00+n+1] |= plane_val;
2056 b[s00+n ] |= plane_val;
2057 b[s00 ] |= plane_val;
2058 break;
2059 case(12):
2060 b[s00+1] |= plane_val;
2061 b[s00 ] |= plane_val;
2062 break;
2063 case(13):
2064 b[s00+n+1] |= plane_val;
2065 b[s00+1] |= plane_val;
2066 b[s00 ] |= plane_val;
2067 break;
2068 case(14):
2069 b[s00+n ] |= plane_val;
2070 b[s00+1] |= plane_val;
2071 b[s00 ] |= plane_val;
2072 break;
2073 case(15):
2074 b[s00+n+1] |= plane_val;
2075 b[s00+n ] |= plane_val;
2076 b[s00+1] |= plane_val;
2077 b[s00 ] |= plane_val;
2078 break;
2079 }
2080
2081 /*
2082 b[s10+1] |= ((LONGLONG) ( a[k] & 1)) << bit;
2083 b[s10 ] |= ((((LONGLONG)a[k])>>1) & 1) << bit;
2084 b[s00+1] |= ((((LONGLONG)a[k])>>2) & 1) << bit;
2085 b[s00 ] |= ((((LONGLONG)a[k])>>3) & 1) << bit;
2086 */
2087 s00 += 2;
2088 /* s10 += 2; */
2089 k += 1;
2090 }
2091 if (j < ny) {
2092 /*
2093 * row size is odd, do last element in row
2094 * s00+1, s10+1 are off edge
2095 */
2096
2097 switch (a[k]) {
2098 case(0):
2099 break;
2100 case(1):
2101 break;
2102 case(2):
2103 b[s00+n ] |= plane_val;
2104 break;
2105 case(3):
2106 b[s00+n ] |= plane_val;
2107 break;
2108 case(4):
2109 break;
2110 case(5):
2111 break;
2112 case(6):
2113 b[s00+n ] |= plane_val;
2114 break;
2115 case(7):
2116 b[s00+n ] |= plane_val;
2117 break;
2118 case(8):
2119 b[s00 ] |= plane_val;
2120 break;
2121 case(9):
2122 b[s00 ] |= plane_val;
2123 break;
2124 case(10):
2125 b[s00+n ] |= plane_val;
2126 b[s00 ] |= plane_val;
2127 break;
2128 case(11):
2129 b[s00+n ] |= plane_val;
2130 b[s00 ] |= plane_val;
2131 break;
2132 case(12):
2133 b[s00 ] |= plane_val;
2134 break;
2135 case(13):
2136 b[s00 ] |= plane_val;
2137 break;
2138 case(14):
2139 b[s00+n ] |= plane_val;
2140 b[s00 ] |= plane_val;
2141 break;
2142 case(15):
2143 b[s00+n ] |= plane_val;
2144 b[s00 ] |= plane_val;
2145 break;
2146 }
2147 /*
2148 b[s10 ] |= ((((LONGLONG)a[k])>>1) & 1) << bit;
2149 b[s00 ] |= ((((LONGLONG)a[k])>>3) & 1) << bit;
2150 */
2151 k += 1;
2152 }
2153 }
2154 if (i < nx) {
2155 /*
2156 * column size is odd, do last row
2157 * s10, s10+1 are off edge
2158 */
2159 s00 = n*i;
2160 for (j = 0; j<ny-1; j += 2) {
2161
2162 switch (a[k]) {
2163 case(0):
2164 break;
2165 case(1):
2166 break;
2167 case(2):
2168 break;
2169 case(3):
2170 break;
2171 case(4):
2172 b[s00+1] |= plane_val;
2173 break;
2174 case(5):
2175 b[s00+1] |= plane_val;
2176 break;
2177 case(6):
2178 b[s00+1] |= plane_val;
2179 break;
2180 case(7):
2181 b[s00+1] |= plane_val;
2182 break;
2183 case(8):
2184 b[s00 ] |= plane_val;
2185 break;
2186 case(9):
2187 b[s00 ] |= plane_val;
2188 break;
2189 case(10):
2190 b[s00 ] |= plane_val;
2191 break;
2192 case(11):
2193 b[s00 ] |= plane_val;
2194 break;
2195 case(12):
2196 b[s00+1] |= plane_val;
2197 b[s00 ] |= plane_val;
2198 break;
2199 case(13):
2200 b[s00+1] |= plane_val;
2201 b[s00 ] |= plane_val;
2202 break;
2203 case(14):
2204 b[s00+1] |= plane_val;
2205 b[s00 ] |= plane_val;
2206 break;
2207 case(15):
2208 b[s00+1] |= plane_val;
2209 b[s00 ] |= plane_val;
2210 break;
2211 }
2212
2213 /*
2214 b[s00+1] |= ((((LONGLONG)a[k])>>2) & 1) << bit;
2215 b[s00 ] |= ((((LONGLONG)a[k])>>3) & 1) << bit;
2216 */
2217 s00 += 2;
2218 k += 1;
2219 }
2220 if (j < ny) {
2221 /*
2222 * both row and column size are odd, do corner element
2223 * s00+1, s10, s10+1 are off edge
2224 */
2225
2226 switch (a[k]) {
2227 case(0):
2228 break;
2229 case(1):
2230 break;
2231 case(2):
2232 break;
2233 case(3):
2234 break;
2235 case(4):
2236 break;
2237 case(5):
2238 break;
2239 case(6):
2240 break;
2241 case(7):
2242 break;
2243 case(8):
2244 b[s00 ] |= plane_val;
2245 break;
2246 case(9):
2247 b[s00 ] |= plane_val;
2248 break;
2249 case(10):
2250 b[s00 ] |= plane_val;
2251 break;
2252 case(11):
2253 b[s00 ] |= plane_val;
2254 break;
2255 case(12):
2256 b[s00 ] |= plane_val;
2257 break;
2258 case(13):
2259 b[s00 ] |= plane_val;
2260 break;
2261 case(14):
2262 b[s00 ] |= plane_val;
2263 break;
2264 case(15):
2265 b[s00 ] |= plane_val;
2266 break;
2267 }
2268 /*
2269 b[s00 ] |= ((((LONGLONG)a[k])>>3) & 1) << bit;
2270 */
2271 k += 1;
2272 }
2273 }
2274 }
2275
2276 /* ############################################################################ */
2277 static void
read_bdirect(unsigned char * infile,int a[],int n,int nqx,int nqy,unsigned char scratch[],int bit)2278 read_bdirect(unsigned char *infile, int a[], int n, int nqx, int nqy, unsigned char scratch[], int bit)
2279 {
2280 /*
2281 * read bit image packed 4 pixels/nybble
2282 */
2283 /*
2284 int i;
2285 for (i = 0; i < ((nqx+1)/2) * ((nqy+1)/2); i++) {
2286 scratch[i] = input_nybble(infile);
2287 }
2288 */
2289 input_nnybble(infile, ((nqx+1)/2) * ((nqy+1)/2), scratch);
2290
2291 /*
2292 * insert in bitplane BIT of image A
2293 */
2294 qtree_bitins(scratch,nqx,nqy,a,n,bit);
2295 }
2296 /* ############################################################################ */
2297 static void
read_bdirect64(unsigned char * infile,LONGLONG a[],int n,int nqx,int nqy,unsigned char scratch[],int bit)2298 read_bdirect64(unsigned char *infile, LONGLONG a[], int n, int nqx, int nqy, unsigned char scratch[], int bit)
2299 {
2300 /*
2301 * read bit image packed 4 pixels/nybble
2302 */
2303 /*
2304 int i;
2305 for (i = 0; i < ((nqx+1)/2) * ((nqy+1)/2); i++) {
2306 scratch[i] = input_nybble(infile);
2307 }
2308 */
2309 input_nnybble(infile, ((nqx+1)/2) * ((nqy+1)/2), scratch);
2310
2311 /*
2312 * insert in bitplane BIT of image A
2313 */
2314 qtree_bitins64(scratch,nqx,nqy,a,n,bit);
2315 }
2316
2317 /* ############################################################################ */
2318 /*
2319 * Huffman decoding for fixed codes
2320 *
2321 * Coded values range from 0-15
2322 *
2323 * Huffman code values (hex):
2324 *
2325 * 3e, 00, 01, 08, 02, 09, 1a, 1b,
2326 * 03, 1c, 0a, 1d, 0b, 1e, 3f, 0c
2327 *
2328 * and number of bits in each code:
2329 *
2330 * 6, 3, 3, 4, 3, 4, 5, 5,
2331 * 3, 5, 4, 5, 4, 5, 6, 4
2332 */
input_huffman(unsigned char * infile)2333 static int input_huffman(unsigned char *infile)
2334 {
2335 int c;
2336
2337 /*
2338 * get first 3 bits to start
2339 */
2340 c = input_nbits(infile,3);
2341 if (c < 4) {
2342 /*
2343 * this is all we need
2344 * return 1,2,4,8 for c=0,1,2,3
2345 */
2346 return(1<<c);
2347 }
2348 /*
2349 * get the next bit
2350 */
2351 c = input_bit(infile) | (c<<1);
2352 if (c < 13) {
2353 /*
2354 * OK, 4 bits is enough
2355 */
2356 switch (c) {
2357 case 8 : return(3);
2358 case 9 : return(5);
2359 case 10 : return(10);
2360 case 11 : return(12);
2361 case 12 : return(15);
2362 }
2363 }
2364 /*
2365 * get yet another bit
2366 */
2367 c = input_bit(infile) | (c<<1);
2368 if (c < 31) {
2369 /*
2370 * OK, 5 bits is enough
2371 */
2372 switch (c) {
2373 case 26 : return(6);
2374 case 27 : return(7);
2375 case 28 : return(9);
2376 case 29 : return(11);
2377 case 30 : return(13);
2378 }
2379 }
2380 /*
2381 * need the 6th bit
2382 */
2383 c = input_bit(infile) | (c<<1);
2384 if (c == 62) {
2385 return(0);
2386 } else {
2387 return(14);
2388 }
2389 }
2390
2391 /* ############################################################################ */
2392 /* ############################################################################ */
2393 /* Copyright (c) 1993 Association of Universities for Research
2394 * in Astronomy. All rights reserved. Produced under National
2395 * Aeronautics and Space Administration Contract No. NAS5-26555.
2396 */
2397 /* qread.c Read binary data
2398 *
2399 * Programmer: R. White Date: 11 March 1991
2400 */
2401
readint(unsigned char * infile)2402 static int readint(unsigned char *infile)
2403 {
2404 int a,i;
2405 unsigned char b[4];
2406
2407 /* Read integer A one byte at a time from infile.
2408 *
2409 * This is portable from Vax to Sun since it eliminates the
2410 * need for byte-swapping.
2411 *
2412 * This routine is only called to read the first 3 values
2413 * in the compressed file, so it doesn't have to be
2414 * super-efficient
2415 */
2416 for (i=0; i<4; i++) qread(infile,(char *) &b[i],1);
2417 a = b[0];
2418 for (i=1; i<4; i++) a = (a<<8) + b[i];
2419 return(a);
2420 }
2421
2422 /* ############################################################################ */
readlonglong(unsigned char * infile)2423 static LONGLONG readlonglong(unsigned char *infile)
2424 {
2425 int i;
2426 LONGLONG a;
2427 unsigned char b[8];
2428
2429 /* Read integer A one byte at a time from infile.
2430 *
2431 * This is portable from Vax to Sun since it eliminates the
2432 * need for byte-swapping.
2433 *
2434 * This routine is only called to read the first 3 values
2435 * in the compressed file, so it doesn't have to be
2436 * super-efficient
2437 */
2438 for (i=0; i<8; i++) qread(infile,(char *) &b[i],1);
2439 a = b[0];
2440 for (i=1; i<8; i++) a = (a<<8) + b[i];
2441 return(a);
2442 }
2443
2444 /* ############################################################################ */
qread(unsigned char * file,char buffer[],int n)2445 static void qread(unsigned char *file, char buffer[], int n)
2446 {
2447 /*
2448 * read n bytes from file into buffer
2449 *
2450 */
2451
2452 memcpy(buffer, &file[nextchar], n);
2453 nextchar += n;
2454 }
2455
2456 /* ############################################################################ */
2457 /* ############################################################################ */
2458 /* Copyright (c) 1993 Association of Universities for Research
2459 * in Astronomy. All rights reserved. Produced under National
2460 * Aeronautics and Space Administration Contract No. NAS5-26555.
2461 */
2462
2463 /* BIT INPUT ROUTINES */
2464
2465 /* THE BIT BUFFER */
2466
2467 static int buffer2; /* Bits waiting to be input */
2468 static int bits_to_go; /* Number of bits still in buffer */
2469
2470 /* INITIALIZE BIT INPUT */
2471
2472 /* ############################################################################ */
start_inputing_bits(void)2473 static void start_inputing_bits(void)
2474 {
2475 /*
2476 * Buffer starts out with no bits in it
2477 */
2478 bits_to_go = 0;
2479 }
2480
2481 /* ############################################################################ */
2482 /* INPUT A BIT */
2483
input_bit(unsigned char * infile)2484 static int input_bit(unsigned char *infile)
2485 {
2486 if (bits_to_go == 0) { /* Read the next byte if no */
2487
2488 buffer2 = infile[nextchar];
2489 nextchar++;
2490
2491 bits_to_go = 8;
2492 }
2493 /*
2494 * Return the next bit
2495 */
2496 bits_to_go -= 1;
2497 return((buffer2>>bits_to_go) & 1);
2498 }
2499
2500 /* ############################################################################ */
2501 /* INPUT N BITS (N must be <= 8) */
2502
input_nbits(unsigned char * infile,int n)2503 static int input_nbits(unsigned char *infile, int n)
2504 {
2505 /* AND mask for retreiving the right-most n bits */
2506 static int mask[9] = {0, 1, 3, 7, 15, 31, 63, 127, 255};
2507
2508 if (bits_to_go < n) {
2509 /*
2510 * need another byte's worth of bits
2511 */
2512
2513 buffer2 = (buffer2<<8) | (int) infile[nextchar];
2514 nextchar++;
2515 bits_to_go += 8;
2516 }
2517 /*
2518 * now pick off the first n bits
2519 */
2520 bits_to_go -= n;
2521
2522 /* there was a slight gain in speed by replacing the following line */
2523 /* return( (buffer2>>bits_to_go) & ((1<<n)-1) ); */
2524 return( (buffer2>>bits_to_go) & (*(mask+n)) );
2525 }
2526 /* ############################################################################ */
2527 /* INPUT 4 BITS */
2528
input_nybble(unsigned char * infile)2529 static int input_nybble(unsigned char *infile)
2530 {
2531 if (bits_to_go < 4) {
2532 /*
2533 * need another byte's worth of bits
2534 */
2535
2536 buffer2 = (buffer2<<8) | (int) infile[nextchar];
2537 nextchar++;
2538 bits_to_go += 8;
2539 }
2540 /*
2541 * now pick off the first 4 bits
2542 */
2543 bits_to_go -= 4;
2544
2545 return( (buffer2>>bits_to_go) & 15 );
2546 }
2547 /* ############################################################################ */
2548 /* INPUT array of 4 BITS */
2549
input_nnybble(unsigned char * infile,int n,unsigned char array[])2550 static int input_nnybble(unsigned char *infile, int n, unsigned char array[])
2551 {
2552 /* copy n 4-bit nybbles from infile to the lower 4 bits of array */
2553
2554 int ii, kk, shift1, shift2;
2555
2556 /* forcing byte alignment doesn;t help, and even makes it go slightly slower
2557 if (bits_to_go != 8) input_nbits(infile, bits_to_go);
2558 */
2559 if (n == 1) {
2560 array[0] = input_nybble(infile);
2561 return(0);
2562 }
2563
2564 if (bits_to_go == 8) {
2565 /*
2566 already have 2 full nybbles in buffer2, so
2567 backspace the infile array to reuse last char
2568 */
2569 nextchar--;
2570 bits_to_go = 0;
2571 }
2572
2573 /* bits_to_go now has a value in the range 0 - 7. After adding */
2574 /* another byte, bits_to_go effectively will be in range 8 - 15 */
2575
2576 shift1 = bits_to_go + 4; /* shift1 will be in range 4 - 11 */
2577 shift2 = bits_to_go; /* shift2 will be in range 0 - 7 */
2578 kk = 0;
2579
2580 /* special case */
2581 if (bits_to_go == 0)
2582 {
2583 for (ii = 0; ii < n/2; ii++) {
2584 /*
2585 * refill the buffer with next byte
2586 */
2587 buffer2 = (buffer2<<8) | (int) infile[nextchar];
2588 nextchar++;
2589 array[kk] = (int) ((buffer2>>4) & 15);
2590 array[kk + 1] = (int) ((buffer2) & 15); /* no shift required */
2591 kk += 2;
2592 }
2593 }
2594 else
2595 {
2596 for (ii = 0; ii < n/2; ii++) {
2597 /*
2598 * refill the buffer with next byte
2599 */
2600 buffer2 = (buffer2<<8) | (int) infile[nextchar];
2601 nextchar++;
2602 array[kk] = (int) ((buffer2>>shift1) & 15);
2603 array[kk + 1] = (int) ((buffer2>>shift2) & 15);
2604 kk += 2;
2605 }
2606 }
2607
2608
2609 if (ii * 2 != n) { /* have to read last odd byte */
2610 array[n-1] = input_nybble(infile);
2611 }
2612
2613 return( (buffer2>>bits_to_go) & 15 );
2614 }
2615