1 package nom.tam.fits.compression.algorithm.hcompress;
2
3 import java.nio.ByteBuffer;
4
5 /*
6 * #%L
7 * nom.tam FITS library
8 * %%
9 * Copyright (C) 1996 - 2024 nom-tam-fits
10 * %%
11 * This is free and unencumbered software released into the public domain.
12 *
13 * Anyone is free to copy, modify, publish, use, compile, sell, or
14 * distribute this software, either in source code form or as a compiled
15 * binary, for any purpose, commercial or non-commercial, and by any
16 * means.
17 *
18 * In jurisdictions that recognize copyright laws, the author or authors
19 * of this software dedicate any and all copyright interest in the
20 * software to the public domain. We make this dedication for the benefit
21 * of the public at large and to the detriment of our heirs and
22 * successors. We intend this dedication to be an overt act of
23 * relinquishment in perpetuity of all present and future rights to this
24 * software under copyright law.
25 *
26 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
27 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
28 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
29 * IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR
30 * OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
31 * ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
32 * OTHER DEALINGS IN THE SOFTWARE.
33 * #L%
34 */
35
36 import static nom.tam.fits.compression.algorithm.hcompress.HCompress.BITS_OF_1_BYTE;
37 import static nom.tam.fits.compression.algorithm.hcompress.HCompress.BITS_OF_1_NYBBLE;
38 import static nom.tam.fits.compression.algorithm.hcompress.HCompress.BYTE_MASK;
39 import static nom.tam.fits.compression.algorithm.hcompress.HCompress.NYBBLE_MASK;
40 import static nom.tam.fits.compression.algorithm.hcompress.HCompress.ROUNDING_HALF;
41
42 /**
43 * (<i>for internal use</i>) A hierarchical data decompression algoritm, e.g. for the Hubble Data Archive and the STScI
44 * Digital Sky Survey.
45 * <p>
46 * The original decompression code was written by R. White at the STScI and included (ported to c and adapted) in
47 * cfitsio by William Pence, NASA/GSFC. That code was then ported to java by R. van Nieuwenhoven. Later it was massively
48 * refactored to harmonize the different compression algorithms and reduce the duplicate code pieces without obscuring
49 * the algorithm itself as far as possible. The original site for the algorithm is
50 * </p>
51 * <p>
52 * See <a href="http://www.stsci.edu/software/hcompress.html">http://www.stsci.edu/software/hcompress.html</a>
53 * </p>
54 *
55 * @author Richard White
56 * @author William Pence
57 * @author Richard van Nieuwenhoven
58 */
59 public class HDecompress {
60
61 private static class LongArrayPointer {
62
63 private final long[] a;
64
65 private int offset;
66
67 LongArrayPointer(long[] tmp) {
68 a = tmp;
69 offset = 0;
70 }
71
72 public void bitOr(int i, long planeVal) {
73 a[offset + i] |= planeVal;
74
75 }
76
77 public LongArrayPointer copy(int extraOffset) {
78 LongArrayPointer intAP = new LongArrayPointer(a);
79 intAP.offset = offset + extraOffset;
80 return intAP;
81 }
82
83 public long get() {
84 return a[offset];
85 }
86
87 public long get(int i) {
88 return a[offset + i];
89 }
90
91 public void set(int i, long value) {
92 a[offset + i] = value;
93
94 }
95
96 public void set(long value) {
97 a[offset] = value;
98 }
99
100 }
101
102 private static final byte[] CODE_MAGIC = {(byte) 0xDD, (byte) 0x99};
103
104 private static final int[] MASKS = {0, 1, 3, 7, 15, 31, 63, 127, 255};
105
106 private static final byte ZERO = 0;
107
108 private static final byte BIT_ONE = 1;
109
110 private static final byte BIT_TWO = 2;
111
112 private static final byte BIT_THREE = 4;
113
114 private static final byte BIT_FOUR = 8;
115
116 /**
117 * these N constants are obscuring the algorithm and should get some explaining javadoc if somebody understands the
118 * algorithm.
119 */
120 private static final int N03 = 3;
121
122 private static final int N04 = 4;
123
124 private static final int N05 = 5;
125
126 private static final int N06 = 6;
127
128 private static final int N07 = 7;
129
130 private static final int N08 = 8;
131
132 private static final int N09 = 9;
133
134 private static final int N10 = 10;
135
136 private static final int N11 = 11;
137
138 private static final int N12 = 12;
139
140 private static final int N13 = 13;
141
142 private static final int N14 = 14;
143
144 private static final int N15 = 15;
145
146 private static final int N26 = 26;
147
148 private static final int N27 = 27;
149
150 private static final int N28 = 28;
151
152 private static final int N29 = 29;
153
154 private static final int N30 = 30;
155
156 private static final int N31 = 31;
157
158 private static final int N62 = 62;
159
160 private static final int N63 = 63;
161
162 /**
163 * Number of bits still in buffer
164 */
165 private int bitsToGo;
166
167 /** Bits waiting to be input */
168 private int buffer2;
169
170 private int nx;
171
172 private int ny;
173
174 private int scale;
175
176 /**
177 * log2n is log2 of max(nx,ny) rounded up to next power of 2
178 */
179 private int calculateLog2N(int nmax) {
180 int log2n;
181 log2n = (int) (Math.log(nmax) / Math.log(2.0) + ROUNDING_HALF);
182 if (nmax > 1 << log2n) {
183 log2n++;
184 }
185 return log2n;
186 }
187
188 /**
189 * char *infile; input file long *a; address of output tiledImageOperation [nx][ny] int *nx,*ny; size of output
190 * tiledImageOperation int *scale; scale factor for digitization
191 *
192 * @param infile
193 * @param a
194 */
195 private void decode64(ByteBuffer infile, LongArrayPointer a) {
196 byte[] nbitplanes = new byte[N03];
197 byte[] tmagic = new byte[2];
198
199 /*
200 * File starts either with special 2-byte magic code or with FITS keyword "SIMPLE ="
201 */
202 infile.get(tmagic);
203 /*
204 * check for correct magic code value
205 */
206 if (tmagic[0] != CODE_MAGIC[0] || tmagic[1] != CODE_MAGIC[1]) {
207 throw new RuntimeException("Compression error");
208 }
209 nx = infile.getInt(); /* x size of image */
210 ny = infile.getInt(); /* y size of image */
211 scale = infile.getInt(); /* scale factor for digitization */
212
213 /* sum of all pixels */
214 long sumall = infile.getLong();
215 /* # bits in quadrants */
216
217 infile.get(nbitplanes);
218
219 dodecode64(infile, a, nbitplanes);
220 /*
221 * put sum of all pixels back into pixel 0
222 */
223 a.set(0, sumall);
224 }
225
226 /**
227 * decompress the input byte stream using the H-compress algorithm input - input tiledImageOperation of compressed
228 * bytes a - pre-allocated tiledImageOperation to hold the output uncompressed image nx - returned X axis size ny -
229 * returned Y axis size NOTE: the nx and ny dimensions as defined within this code are reversed from the usual FITS
230 * notation. ny is the fastest varying dimension, which is usually considered the X axis in the FITS image display
231 *
232 * @param input the input buffer to decompress
233 * @param smooth should the image be smoothed
234 * @param aa the resulting long tiledImageOperation
235 */
236 public void decompress(ByteBuffer input, boolean smooth, long[] aa) {
237
238 LongArrayPointer a = new LongArrayPointer(aa);
239
240 /* decode the input tiledImageOperation */
241
242 decode64(input, a);
243
244 /*
245 * Un-Digitize
246 */
247 undigitize64(a);
248
249 /*
250 * Inverse H-transform
251 */
252 hinv64(a, smooth);
253
254 }
255
256 /**
257 * long a[]; int nx,ny; Array dimensions are [nx][ny] unsigned char nbitplanes[3]; Number of bit planes in quadrants
258 */
259 private int dodecode64(ByteBuffer infile, LongArrayPointer a, byte[] nbitplanes) {
260 int nel = nx * ny;
261 int nx2 = (nx + 1) / 2;
262 int ny2 = (ny + 1) / 2;
263 /*
264 * initialize a to zero
265 */
266 for (int i = 0; i < nel; i++) {
267 a.set(i, 0);
268 }
269 /*
270 * Initialize bit input
271 */
272 startInputingBits();
273 /*
274 * read bit planes for each quadrant
275 */
276 qtreeDecode64(infile, a.copy(0), ny, nx2, ny2, nbitplanes[0]);
277
278 qtreeDecode64(infile, a.copy(ny2), ny, nx2, ny / 2, nbitplanes[1]);
279
280 qtreeDecode64(infile, a.copy(ny * nx2), ny, nx / 2, ny2, nbitplanes[1]);
281
282 qtreeDecode64(infile, a.copy(ny * nx2 + ny2), ny, nx / 2, ny / 2, nbitplanes[2]);
283
284 /*
285 * make sure there is an EOF symbol (nybble=0) at end
286 */
287 if (inputNybble(infile) != 0) {
288 throw new RuntimeException("Compression error");
289 }
290 /*
291 * now get the sign bits Re-initialize bit input
292 */
293 startInputingBits();
294 for (int i = 0; i < nel; i++) {
295 if (a.get(i) != 0) {
296 if (inputBit(infile) != 0) {
297 a.set(i, -a.get(i));
298 }
299 }
300 }
301 return 0;
302 }
303
304 /**
305 * int smooth; 0 for no smoothing, else smooth during inversion int scale; used if smoothing is specified
306 */
307 private int hinv64(LongArrayPointer a, boolean smooth) {
308 int nmax = nx > ny ? nx : ny;
309 int log2n = calculateLog2N(nmax);
310 // get temporary storage for shuffling elements
311 long[] tmp = new long[(nmax + 1) / 2];
312 // set up masks, rounding parameters
313 int shift = 1;
314 long bit0 = (long) 1 << log2n - 1;
315 long bit1 = bit0 << 1;
316 long bit2 = bit0 << 2;
317 long mask0 = -bit0;
318 long mask1 = mask0 << 1;
319 long mask2 = mask0 << 2;
320 long prnd0 = bit0 >> 1;
321 long prnd1 = bit1 >> 1;
322 long prnd2 = bit2 >> 1;
323 long nrnd0 = prnd0 - 1;
324 long nrnd1 = prnd1 - 1;
325 long nrnd2 = prnd2 - 1;
326 // round h0 to multiple of bit2
327 a.set(0, a.get(0) + (a.get(0) >= 0 ? prnd2 : nrnd2) & mask2);
328 // do log2n expansions We're indexing a as a 2-D tiledImageOperation
329 // with dimensions
330 // (nx,ny).
331 int nxtop = 1;
332 int nytop = 1;
333 int nxf = nx;
334 int nyf = ny;
335 int c = 1 << log2n;
336 int i;
337 for (int k = log2n - 1; k >= 0; k--) {
338 // this somewhat cryptic code generates the sequence ntop[k-1] =
339 // (ntop[k]+1)/2, where ntop[log2n] = n
340 c = c >> 1;
341 nxtop = nxtop << 1;
342 nytop = nytop << 1;
343 if (nxf <= c) {
344 nxtop--;
345 } else {
346 nxf -= c;
347 }
348 if (nyf <= c) {
349 nytop--;
350 } else {
351 nyf -= c;
352 }
353 // double shift and fix nrnd0 (because prnd0=0) on last pass
354 if (k == 0) {
355 nrnd0 = 0;
356 shift = 2;
357 }
358 // unshuffle in each dimension to interleave coefficients
359 for (i = 0; i < nxtop; i++) {
360 unshuffle64(a.copy(ny * i), nytop, 1, tmp);
361 }
362 for (int j = 0; j < nytop; j++) {
363 unshuffle64(a.copy(j), nxtop, ny, tmp);
364 }
365 // smooth by interpolating coefficients if SMOOTH != 0
366 if (smooth) {
367 hsmooth64(a, nxtop, nytop);
368 }
369 int oddx = nxtop % 2;
370 int oddy = nytop % 2;
371 for (i = 0; i < nxtop - oddx; i += 2) {
372 int s00 = ny * i; /* s00 is index of a[i,j] */
373 int s10 = s00 + ny; /* s10 is index of a[i+1,j] */
374 for (int j = 0; j < nytop - oddy; j += 2) {
375 long h0 = a.get(s00);
376 long hx = a.get(s10);
377 long hy = a.get(s00 + 1);
378 long hc = a.get(s10 + 1);
379 // round hx and hy to multiple of bit1, hc to multiple of
380 // bit0 h0 is already a multiple of bit2
381 hx = hx + (hx >= 0 ? prnd1 : nrnd1) & mask1;
382 hy = hy + (hy >= 0 ? prnd1 : nrnd1) & mask1;
383 hc = hc + (hc >= 0 ? prnd0 : nrnd0) & mask0;
384 // propagate bit0 of hc to hx,hy
385 long lowbit0 = hc & bit0;
386 hx = hx >= 0 ? hx - lowbit0 : hx + lowbit0;
387 hy = hy >= 0 ? hy - lowbit0 : hy + lowbit0;
388 // Propagate bits 0 and 1 of hc,hx,hy to h0. This could be
389 // simplified if we assume h0>0, but then the inversion
390 // would not be lossless for images with negative pixels.
391 long lowbit1 = (hc ^ hx ^ hy) & bit1;
392 h0 = h0 >= 0 ? h0 + lowbit0 - lowbit1 : h0 + (lowbit0 == 0 ? lowbit1 : lowbit0 - lowbit1);
393 // Divide sums by 2 (4 last time)
394 a.set(s10 + 1, h0 + hx + hy + hc >> shift);
395 a.set(s10, h0 + hx - hy - hc >> shift);
396 a.set(s00 + 1, h0 - hx + hy - hc >> shift);
397 a.set(s00, h0 - hx - hy + hc >> shift);
398 s00 += 2;
399 s10 += 2;
400 }
401 if (oddy != 0) {
402 // do last element in row if row length is odd s00+1, s10+1
403 // are off edge
404 long h0 = a.get(s00);
405 long hx = a.get(s10);
406 hx = (hx >= 0 ? hx + prnd1 : hx + nrnd1) & mask1;
407 long lowbit1 = hx & bit1;
408 h0 = h0 >= 0 ? h0 - lowbit1 : h0 + lowbit1;
409 a.set(s10, h0 + hx >> shift);
410 a.set(s00, h0 - hx >> shift);
411 }
412 }
413 if (oddx != 0) {
414 // do last row if column length is odd s10, s10+1 are off edge
415 int s00 = ny * i;
416 for (int j = 0; j < nytop - oddy; j += 2) {
417 long h0 = a.get(s00);
418 long hy = a.get(s00 + 1);
419 hy = (hy >= 0 ? hy + prnd1 : hy + nrnd1) & mask1;
420 long lowbit1 = hy & bit1;
421 h0 = h0 >= 0 ? h0 - lowbit1 : h0 + lowbit1;
422 a.set(s00 + 1, h0 + hy >> shift);
423 a.set(s00, h0 - hy >> shift);
424 s00 += 2;
425 }
426 if (oddy != 0) {
427 // do corner element if both row and column lengths are odd
428 // s00+1, s10, s10+1 are off edge
429 long h0 = a.get(s00);
430 a.set(s00, h0 >> shift);
431 }
432 }
433 // divide all the masks and rounding values by 2
434 bit1 = bit0;
435 bit0 = bit0 >> 1;
436 mask1 = mask0;
437 mask0 = mask0 >> 1;
438 prnd1 = prnd0;
439 prnd0 = prnd0 >> 1;
440 nrnd1 = nrnd0;
441 nrnd0 = prnd0 - 1;
442 }
443 return 0;
444 }
445
446 /**
447 * long a[]; tiledImageOperation of H-transform coefficients int nxtop,nytop; size of coefficient block to use int
448 * ny; actual 1st dimension of tiledImageOperation int scale; truncation scale factor that was used
449 */
450 private void hsmooth64(LongArrayPointer a, int nxtop, int nytop) {
451 int i, j;
452 int ny2, s10, s00;
453 long hm, h0, hp, hmm, hpm, hmp, hpp, hx2, hy2, diff, dmax, dmin, s, smax, m1, m2;
454
455 /*
456 * Maximum change in coefficients is determined by scale factor. Since we rounded during division (see
457 * digitize.c), the biggest permitted change is scale/2.
458 */
459 smax = scale >> 1;
460 if (smax <= 0) {
461 return;
462 }
463 ny2 = ny << 1;
464 /*
465 * We're indexing a as a 2-D tiledImageOperation with dimensions (nxtop,ny) of which only (nxtop,nytop) are
466 * used. The coefficients on the edge of the tiledImageOperation are not adjusted (which is why the loops below
467 * start at 2 instead of 0 and end at nxtop-2 instead of nxtop.)
468 */
469 /*
470 * Adjust x difference hx
471 */
472 for (i = 2; i < nxtop - 2; i += 2) {
473 s00 = ny * i; /* s00 is index of a[i,j] */
474 s10 = s00 + ny; /* s10 is index of a[i+1,j] */
475 for (j = 0; j < nytop; j += 2) {
476 /*
477 * hp is h0 (mean value) in next x zone, hm is h0 in previous x zone
478 */
479 hm = a.get(s00 - ny2);
480 h0 = a.get(s00);
481 hp = a.get(s00 + ny2);
482 /*
483 * diff = 8 * hx slope that would match h0 in neighboring zones
484 */
485 diff = hp - hm;
486 /*
487 * monotonicity constraints on diff
488 */
489 dmax = Math.max(Math.min(hp - h0, h0 - hm), 0) << 2;
490 dmin = Math.min(Math.max(hp - h0, h0 - hm), 0) << 2;
491 /*
492 * if monotonicity would set slope = 0 then don't change hx. note dmax>=0, dmin<=0.
493 */
494 if (dmin < dmax) {
495 diff = Math.max(Math.min(diff, dmax), dmin);
496 /*
497 * Compute change in slope limited to range +/- smax. Careful with rounding negative numbers when
498 * using shift for divide by 8.
499 */
500 s = diff - (a.get(s10) << N03);
501 s = s >= 0 ? s >> N03 : s + N07 >> N03;
502 s = Math.max(Math.min(s, smax), -smax);
503 a.set(s10, a.get(s10) + s);
504 }
505 s00 += 2;
506 s10 += 2;
507 }
508 }
509 /*
510 * Adjust y difference hy
511 */
512 for (i = 0; i < nxtop; i += 2) {
513 s00 = ny * i + 2;
514 s10 = s00 + ny;
515 for (j = 2; j < nytop - 2; j += 2) {
516 hm = a.get(s00 - 2);
517 h0 = a.get(s00);
518 hp = a.get(s00 + 2);
519 diff = hp - hm;
520 dmax = Math.max(Math.min(hp - h0, h0 - hm), 0) << 2;
521 dmin = Math.min(Math.max(hp - h0, h0 - hm), 0) << 2;
522 if (dmin < dmax) {
523 diff = Math.max(Math.min(diff, dmax), dmin);
524 s = diff - (a.get(s00 + 1) << N03);
525 s = s >= 0 ? s >> N03 : s + N07 >> N03;
526 s = Math.max(Math.min(s, smax), -smax);
527 a.set(s00 + 1, a.get(s00 + 1) + s);
528 }
529 s00 += 2;
530 s10 += 2;
531 }
532 }
533 /*
534 * Adjust curvature difference hc
535 */
536 for (i = 2; i < nxtop - 2; i += 2) {
537 s00 = ny * i + 2;
538 s10 = s00 + ny;
539 for (j = 2; j < nytop - 2; j += 2) {
540 /*
541 * ------------------ y | hmp | | hpp | | ------------------ | | | h0 | | | ------------------ -------x
542 * | hmm | | hpm | ------------------
543 */
544 hmm = a.get(s00 - ny2 - 2);
545 hpm = a.get(s00 + ny2 - 2);
546 hmp = a.get(s00 - ny2 + 2);
547 hpp = a.get(s00 + ny2 + 2);
548 h0 = a.get(s00);
549 /*
550 * diff = 64 * hc value that would match h0 in neighboring zones
551 */
552 diff = hpp + hmm - hmp - hpm;
553 /*
554 * 2 times x,y slopes in this zone
555 */
556 hx2 = a.get(s10) << 1;
557 hy2 = a.get(s00 + 1) << 1;
558 /*
559 * monotonicity constraints on diff
560 */
561 m1 = Math.min(Math.max(hpp - h0, 0) - hx2 - hy2, Math.max(h0 - hpm, 0) + hx2 - hy2);
562 m2 = Math.min(Math.max(h0 - hmp, 0) - hx2 + hy2, Math.max(hmm - h0, 0) + hx2 + hy2);
563 dmax = Math.min(m1, m2) << BITS_OF_1_NYBBLE;
564 m1 = Math.max(Math.min(hpp - h0, 0) - hx2 - hy2, Math.min(h0 - hpm, 0) + hx2 - hy2);
565 m2 = Math.max(Math.min(h0 - hmp, 0) - hx2 + hy2, Math.min(hmm - h0, 0) + hx2 + hy2);
566 dmin = Math.max(m1, m2) << BITS_OF_1_NYBBLE;
567 /*
568 * if monotonicity would set slope = 0 then don't change hc. note dmax>=0, dmin<=0.
569 */
570 if (dmin < dmax) {
571 diff = Math.max(Math.min(diff, dmax), dmin);
572 /*
573 * Compute change in slope limited to range +/- smax. Careful with rounding negative numbers when
574 * using shift for divide by 64.
575 */
576 s = diff - (a.get(s10 + 1) << N06);
577 s = s >= 0 ? s >> N06 : s + N63 >> N06;
578 s = Math.max(Math.min(s, smax), -smax);
579 a.set(s10 + 1, a.get(s10 + 1) + s);
580 }
581 s00 += 2;
582 s10 += 2;
583 }
584 }
585 }
586
587 private int inputBit(ByteBuffer infile) {
588 if (bitsToGo == 0) { /* Read the next byte if no */
589
590 buffer2 = infile.get() & BYTE_MASK;
591
592 bitsToGo = BITS_OF_1_BYTE;
593 }
594 /*
595 * Return the next bit
596 */
597 bitsToGo--;
598 return buffer2 >> bitsToGo & 1;
599 }
600
601 /*
602 * Huffman decoding for fixed codes Coded values range from 0-15 Huffman code values (hex): 3e, 00, 01, 08, 02, 09,
603 * 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,
604 * 6, 4
605 */
606 private int inputHuffman(ByteBuffer infile) {
607 int c;
608
609 /*
610 * get first 3 bits to start
611 */
612 c = inputNbits(infile, N03);
613 if (c < N04) {
614 /*
615 * this is all we need return 1,2,4,8 for c=0,1,2,3
616 */
617 return 1 << c;
618 }
619 /*
620 * get the next bit
621 */
622 c = inputBit(infile) | c << 1;
623 if (c < N13) {
624 /*
625 * OK, 4 bits is enough
626 */
627 switch (c) {
628 case N08:
629 return N03;
630 case N09:
631 return N05;
632 case N10:
633 return N10;
634 case N11:
635 return N12;
636 case N12:
637 return N15;
638 default:
639 }
640 }
641 /*
642 * get yet another bit
643 */
644 c = inputBit(infile) | c << 1;
645 if (c < N31) {
646 /*
647 * OK, 5 bits is enough
648 */
649 switch (c) {
650 case N26:
651 return N06;
652 case N27:
653 return N07;
654 case N28:
655 return N09;
656 case N29:
657 return N11;
658 case N30:
659 return N13;
660 default:
661 }
662 }
663 /*
664 * need the 6th bit
665 */
666 c = inputBit(infile) | c << 1;
667 if (c == N62) {
668 return 0;
669 }
670 return N14;
671 }
672
673 private int inputNbits(ByteBuffer infile, int n) {
674 if (bitsToGo < n) {
675 /*
676 * need another byte's worth of bits
677 */
678
679 buffer2 = buffer2 << BITS_OF_1_BYTE | infile.get() & BYTE_MASK;
680 bitsToGo += BITS_OF_1_BYTE;
681 }
682 /*
683 * now pick off the first n bits
684 */
685 bitsToGo -= n;
686
687 /* there was a slight gain in speed by replacing the following line */
688 /* return( (buffer2>>bits_to_go) & ((1<<n)-1) ); */
689 return buffer2 >> bitsToGo & MASKS[n];
690 }
691
692 /* INITIALIZE BIT INPUT */
693
694 private int inputNnybble(ByteBuffer infile, int n, byte[] array) {
695 /*
696 * copy n 4-bit nybbles from infile to the lower 4 bits of tiledImageOperation
697 */
698
699 int ii, kk, shift1, shift2;
700
701 /*
702 * forcing byte alignment doesn;t help, and even makes it go slightly slower if (bits_to_go != 8)
703 * input_nbits(infile, bits_to_go);
704 */
705 if (n == 1) {
706 array[0] = (byte) inputNybble(infile);
707 return 0;
708 }
709
710 if (bitsToGo == BITS_OF_1_BYTE) {
711 /*
712 * already have 2 full nybbles in buffer2, so backspace the infile tiledImageOperation to reuse last char
713 */
714 infile.position(infile.position() - 1);
715 bitsToGo = 0;
716 }
717
718 /* bits_to_go now has a value in the range 0 - 7. After adding */
719 /* another byte, bits_to_go effectively will be in range 8 - 15 */
720
721 shift1 = bitsToGo + BITS_OF_1_NYBBLE; /*
722 * shift1 will be in range 4 - 11
723 */
724 shift2 = bitsToGo; /* shift2 will be in range 0 - 7 */
725 kk = 0;
726
727 /* special case */
728 if (bitsToGo == 0) {
729 for (ii = 0; ii < n / 2; ii++) {
730 /*
731 * refill the buffer with next byte
732 */
733 buffer2 = buffer2 << BITS_OF_1_BYTE | infile.get() & BYTE_MASK;
734 array[kk] = (byte) (buffer2 >> BITS_OF_1_NYBBLE & NYBBLE_MASK);
735 array[kk + 1] = (byte) (buffer2 & NYBBLE_MASK); /*
736 * no shift required
737 */
738 kk += 2;
739 }
740 } else {
741 for (ii = 0; ii < n / 2; ii++) {
742 /*
743 * refill the buffer with next byte
744 */
745 buffer2 = buffer2 << BITS_OF_1_BYTE | infile.get() & BYTE_MASK;
746 array[kk] = (byte) (buffer2 >> shift1 & NYBBLE_MASK);
747 array[kk + 1] = (byte) (buffer2 >> shift2 & NYBBLE_MASK);
748 kk += 2;
749 }
750 }
751
752 if (ii * 2 != n) { /* have to read last odd byte */
753 array[n - 1] = (byte) inputNybble(infile);
754 }
755
756 return buffer2 >> bitsToGo & NYBBLE_MASK;
757 }
758
759 private int inputNybble(ByteBuffer infile) {
760 if (bitsToGo < BITS_OF_1_NYBBLE) {
761 /*
762 * need another byte's worth of bits
763 */
764
765 buffer2 = buffer2 << BITS_OF_1_BYTE | infile.get() & BYTE_MASK;
766 bitsToGo += BITS_OF_1_BYTE;
767 }
768 /*
769 * now pick off the first 4 bits
770 */
771 bitsToGo -= BITS_OF_1_NYBBLE;
772
773 return buffer2 >> bitsToGo & NYBBLE_MASK;
774 }
775
776 /**
777 * 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
778 * bitplane BIT of B. A,B may NOT be same tiledImageOperation (it wouldn't make sense to be inserting bits into the
779 * same tiledImageOperation anyway.)
780 */
781 private void qtreeBitins64(byte[] a, int lnx, int lny, LongArrayPointer b, int n, int bit) {
782 int i, j, s00;
783 long planeVal = 1L << bit;
784 // expand each 2x2 block
785 ByteBuffer k = ByteBuffer.wrap(a); /* k is index of a[i/2,j/2] */
786 for (i = 0; i < lnx - 1; i += 2) {
787 s00 = n * i; /* s00 is index of b[i,j] */
788 // Note: this code appears to run very slightly faster on a 32-bit
789 // linux machine using s00+n rather than the s10 intermediate
790 // variable
791 // s10 = s00+n; *//* s10 is index of b[i+1,j]
792 for (j = 0; j < lny - 1; j += 2) {
793 byte value = k.get();
794 if ((value & BIT_ONE) != ZERO) {
795 b.bitOr(s00 + n + 1, planeVal);
796 }
797 if ((value & BIT_TWO) != ZERO) {
798 b.bitOr(s00 + n, planeVal);
799 }
800 if ((value & BIT_THREE) != ZERO) {
801 b.bitOr(s00 + 1, planeVal);
802 }
803 if ((value & BIT_FOUR) != ZERO) {
804 b.bitOr(s00, planeVal);
805 }
806 // b.bitOr(s10+1, ((LONGLONG) ( a[k] & 1)) << bit; b.bitOr(s10 ,
807 // ((((LONGLONG)a[k])>>1) & 1) << bit; b.bitOr(s00+1,
808 // ((((LONGLONG)a[k])>>2) & 1) << bit; b.bitOr(s00
809 // ,((((LONGLONG)a[k])>>3) & 1) << bit;
810 s00 += 2;
811 /* s10 += 2; */
812 }
813 if (j < lny) {
814 // row size is odd, do last element in row s00+1, s10+1 are off
815 // edge
816 byte value = k.get();
817 if ((value & BIT_TWO) != ZERO) {
818 b.bitOr(s00 + n, planeVal);
819 }
820 if ((value & BIT_FOUR) != ZERO) {
821 b.bitOr(s00, planeVal);
822 }
823 // b.bitOr(s10 , ((((LONGLONG)a[k])>>1) & 1) << bit; b.bitOr(s00
824 // , ((((LONGLONG)a[k])>>3) & 1) << bit;
825 }
826 }
827 if (i < lnx) {
828 // column size is odd, do last row s10, s10+1 are off edge
829 s00 = n * i;
830 for (j = 0; j < lny - 1; j += 2) {
831 byte value = k.get();
832 if ((value & BIT_THREE) != ZERO) {
833 b.bitOr(s00 + 1, planeVal);
834 }
835 if ((value & BIT_FOUR) != ZERO) {
836 b.bitOr(s00, planeVal);
837 } // b.bitOr(s00+1, ((((LONGLONG)a[k])>>2) & 1) << bit;
838 // b.bitOr(s00 , ((((LONGLONG)a[k])>>3) & 1) << bit;
839 s00 += 2;
840 }
841 if (j < lny) {
842 // both row and column size are odd, do corner element s00+1,
843 // s10, s10+1 are off edge
844 if ((k.get() & BIT_FOUR) != ZERO) {
845 b.bitOr(s00, planeVal);
846 }
847 // b.bitOr(s00 , ((((LONGLONG)a[k])>>3) & 1) << bit;
848 }
849 }
850 }
851
852 /**
853 * 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
854 * tiledImageOperation
855 */
856 private void qtreeCopy(byte[] a, int lnx, int lny, byte[] b, int n) {
857 int i, j, k, nx2, ny2;
858 int s00, s10;
859 // first copy 4-bit values to b start at end in case a,b are same
860 // tiledImageOperation
861 nx2 = (lnx + 1) / 2;
862 ny2 = (lny + 1) / 2;
863 k = ny2 * (nx2 - 1) + ny2 - 1; /* k is index of a[i,j] */
864 for (i = nx2 - 1; i >= 0; i--) {
865 s00 = 2 * (n * i + ny2 - 1); /* s00 is index of b[2*i,2*j] */
866 for (j = ny2 - 1; j >= 0; j--) {
867 b[s00] = a[k];
868 k--;
869 s00 -= 2;
870 }
871 }
872 for (i = 0; i < lnx - 1; i += 2) { // now expand each 2x2 block
873 // Note: Unlike the case in qtree_bitins, this code runs faster on a
874 // 32-bit linux machine using the s10 intermediate variable, rather
875 // that using s00+n. Go figure!
876 s00 = n * i; // s00 is index of b[i,j]
877 s10 = s00 + n; // s10 is index of b[i+1,j]
878 for (j = 0; j < lny - 1; j += 2) {
879 b[s10 + 1] = (b[s00] & BIT_ONE) == ZERO ? ZERO : BIT_ONE;
880 b[s10] = (b[s00] & BIT_TWO) == ZERO ? ZERO : BIT_ONE;
881 b[s00 + 1] = (b[s00] & BIT_THREE) == ZERO ? ZERO : BIT_ONE;
882 b[s00] = (b[s00] & BIT_FOUR) == ZERO ? ZERO : BIT_ONE;
883 s00 += 2;
884 s10 += 2;
885 }
886 if (j < lny) {
887 // row size is odd, do last element in row s00+1, s10+1 are off
888 // edge not worth converting this to use 16 case statements
889 b[s10] = (byte) (b[s00] >> 1 & 1);
890 b[s00] = (byte) (b[s00] >> N03 & 1);
891 }
892 }
893 if (i < lnx) {
894 // column size is odd, do last row s10, s10+1 are off edge
895 s00 = n * i;
896 for (j = 0; j < lny - 1; j += 2) {
897 // not worth converting this to use 16 case statements
898 b[s00 + 1] = (byte) (b[s00] >> 2 & 1);
899 b[s00] = (byte) (b[s00] >> N03 & 1);
900 s00 += 2;
901 }
902 if (j < lny) {
903 // both row and column size are odd, do corner element s00+1,
904 // s10, s10+1 are off edge not worth converting this to use 16
905 // case statements
906 b[s00] = (byte) (b[s00] >> N03 & 1);
907 }
908 }
909 }
910
911 /**
912 * char *infile; long a[]; a is 2-D tiledImageOperation with dimensions (n,n) int n; length of full row in a int
913 * nqx; partial length of row to decode int nqy; partial length of column (<=n) int nbitplanes; number of bitplanes
914 * to decode
915 */
916 private int qtreeDecode64(ByteBuffer infile, LongArrayPointer a, int n, int nqx, int nqy, int nbitplanes) {
917 int k, bit, b;
918 int nx2, ny2, nfx, nfy, c;
919 byte[] scratch;
920
921 /*
922 * log2n is log2 of max(nqx,nqy) rounded up to next power of 2
923 */
924 int nqmax = nqx > nqy ? nqx : nqy;
925 int log2n = calculateLog2N(nqmax);
926 /*
927 * allocate scratch tiledImageOperation for working space
928 */
929 int nqx2 = (nqx + 1) / 2;
930 int nqy2 = (nqy + 1) / 2;
931 scratch = new byte[nqx2 * nqy2];
932
933 /*
934 * now decode each bit plane, starting at the top A is assumed to be initialized to zero
935 */
936 for (bit = nbitplanes - 1; bit >= 0; bit--) {
937 /*
938 * Was bitplane was quadtree-coded or written directly?
939 */
940 b = inputNybble(infile);
941
942 if (b == 0) {
943 /*
944 * bit map was written directly
945 */
946 readBdirect64(infile, a, n, nqx, nqy, scratch, bit);
947 } else if (b != NYBBLE_MASK) {
948 throw new RuntimeException("Compression error");
949 } else {
950 /*
951 * bitmap was quadtree-coded, do log2n expansions read first code
952 */
953 scratch[0] = (byte) inputHuffman(infile);
954 /*
955 * now do log2n expansions, reading codes from file as necessary
956 */
957 nx2 = 1;
958 ny2 = 1;
959 nfx = nqx;
960 nfy = nqy;
961 c = 1 << log2n;
962 for (k = 1; k < log2n; k++) {
963 /*
964 * this somewhat cryptic code generates the sequence n[k-1] = (n[k]+1)/2 where n[log2n]=nqx or nqy
965 */
966 c = c >> 1;
967 nx2 = nx2 << 1;
968 ny2 = ny2 << 1;
969 if (nfx <= c) {
970 nx2--;
971 } else {
972 nfx -= c;
973 }
974 if (nfy <= c) {
975 ny2--;
976 } else {
977 nfy -= c;
978 }
979 qtreeExpand(infile, scratch, nx2, ny2, scratch);
980 }
981 /*
982 * now copy last set of 4-bit codes to bitplane bit of tiledImageOperation a
983 */
984 qtreeBitins64(scratch, nqx, nqy, a, n, bit);
985 }
986 }
987 return 0;
988 }
989
990 /*
991 * do one quadtree expansion step on tiledImageOperation a[(nqx+1)/2,(nqy+1)/2] results put into b[nqx,nqy] (which
992 * may be the same as a)
993 */
994 private void qtreeExpand(ByteBuffer infile, byte[] a, int nx2, int ny2, byte[] b) {
995 int i;
996
997 /*
998 * first copy a to b, expanding each 4-bit value
999 */
1000 qtreeCopy(a, nx2, ny2, b, ny2);
1001 /*
1002 * now read new 4-bit values into b for each non-zero element
1003 */
1004 for (i = nx2 * ny2 - 1; i >= 0; i--) {
1005 if (b[i] != 0) {
1006 b[i] = (byte) inputHuffman(infile);
1007 }
1008 }
1009 }
1010
1011 private void readBdirect64(ByteBuffer infile, LongArrayPointer a, int n, int nqx, int nqy, byte[] scratch, int bit) {
1012 /*
1013 * read bit image packed 4 pixels/nybble
1014 */
1015 /*
1016 * int i; for (i = 0; i < ((nqx+1)/2) * ((nqy+1)/2); i++) { scratch[i] = input_nybble(infile); }
1017 */
1018 inputNnybble(infile, (nqx + 1) / 2 * ((nqy + 1) / 2), scratch);
1019
1020 /*
1021 * insert in bitplane BIT of image A
1022 */
1023 qtreeBitins64(scratch, nqx, nqy, a, n, bit);
1024 }
1025
1026 /*
1027 * ########################################################################## ##
1028 */
1029 private void startInputingBits() {
1030 /*
1031 * Buffer starts out with no bits in it
1032 */
1033 bitsToGo = 0;
1034 }
1035
1036 private void undigitize64(LongArrayPointer a) {
1037 long scale64;
1038
1039 /*
1040 * multiply by scale
1041 */
1042 if (scale <= 1) {
1043 return;
1044 }
1045 scale64 = scale; /*
1046 * use a 64-bit int for efficiency in the big loop
1047 */
1048
1049 for (int index = 0; index < a.a.length; index++) {
1050 a.a[index] = a.a[index] * scale64;
1051 }
1052 }
1053
1054 /**
1055 * long a[]; tiledImageOperation to shuffle int n; number of elements to shuffle int n2; second dimension long
1056 * tmp[]; scratch storage
1057 */
1058 private void unshuffle64(LongArrayPointer a, int n, int n2, long[] tmp) {
1059 int i;
1060 int nhalf;
1061 LongArrayPointer p1, p2, pt;
1062
1063 /*
1064 * copy 2nd half of tiledImageOperation to tmp
1065 */
1066 nhalf = n + 1 >> 1;
1067 pt = new LongArrayPointer(tmp);
1068 p1 = a.copy(n2 * nhalf); /* pointer to a[i] */
1069 for (i = nhalf; i < n; i++) {
1070 pt.set(p1.get());
1071 p1.offset += n2;
1072 pt.offset++;
1073 }
1074 /*
1075 * distribute 1st half of tiledImageOperation to even elements
1076 */
1077 p2 = a.copy(n2 * (nhalf - 1)); /* pointer to a[i] */
1078 p1 = a.copy(n2 * (nhalf - 1) << 1); /* pointer to a[2*i] */
1079 for (i = nhalf - 1; i >= 0; i--) {
1080 p1.set(p2.get());
1081 p2.offset -= n2;
1082 p1.offset -= n2 + n2;
1083 }
1084 /*
1085 * now distribute 2nd half of tiledImageOperation (in tmp) to odd elements
1086 */
1087 pt = new LongArrayPointer(tmp);
1088 p1 = a.copy(n2); /* pointer to a[i] */
1089 for (i = 1; i < n; i += 2) {
1090 p1.set(pt.get());
1091 p1.offset += n2 + n2;
1092 pt.offset++;
1093 }
1094 }
1095 }