1 /*
2 The following code was written by Richard White at STScI and made
3 available for use in CFITSIO in July 1999. These routines were
4 originally contained in 2 source files: rcomp.c and rdecomp.c,
5 and the 'include' file now called ricecomp.h was originally called buffer.h.
6
7 Note that beginning with CFITSIO v3.08, EOB checking was removed to improve
8 speed, and so now the input compressed bytes buffers must have been
9 allocated big enough so that they will never be overflowed. A simple
10 rule of thumb that guarantees the buffer will be large enough is to make
11 it 1% larger than the size of the input array of pixels that are being
12 compressed.
13
14 */
15
16 /*----------------------------------------------------------*/
17 /* */
18 /* START OF SOURCE FILE ORIGINALLY CALLED rcomp.c */
19 /* */
20 /*----------------------------------------------------------*/
21 /* @(#) rcomp.c 1.5 99/03/01 12:40:27 */
22 /* rcomp.c Compress image line using
23 * (1) Difference of adjacent pixels
24 * (2) Rice algorithm coding
25 *
26 * Returns number of bytes written to code buffer or
27 * -1 on failure
28 */
29
30 #include <stdio.h>
31 #include <stdlib.h>
32 #include <string.h>
33
34 /*
35 * nonzero_count is lookup table giving number of bits in 8-bit values not including
36 * leading zeros used in fits_rdecomp, fits_rdecomp_short and fits_rdecomp_byte
37 */
38 static const int nonzero_count[256] = {
39 0,
40 1,
41 2, 2,
42 3, 3, 3, 3,
43 4, 4, 4, 4, 4, 4, 4, 4,
44 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
45 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
46 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
47 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
48 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
49 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
50 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
51 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
52 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
53 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
54 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
55 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
56 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
57 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
58 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8};
59
60 typedef unsigned char Buffer_t;
61
62 typedef struct {
63 int bitbuffer; /* bit buffer */
64 int bits_to_go; /* bits to go in buffer */
65 Buffer_t *start; /* start of buffer */
66 Buffer_t *current; /* current position in buffer */
67 Buffer_t *end; /* end of buffer */
68 } Buffer;
69
70 #define putcbuf(c,mf) ((*(mf->current)++ = c), 0)
71
72 #include "fitsio2.h"
73
74 static void start_outputing_bits(Buffer *buffer);
75 static int done_outputing_bits(Buffer *buffer);
76 static int output_nbits(Buffer *buffer, int bits, int n);
77
78 /* only used for diagnoistics
79 static int case1, case2, case3;
80 int fits_get_case(int *c1, int*c2, int*c3) {
81
82 *c1 = case1;
83 *c2 = case2;
84 *c3 = case3;
85 return(0);
86 }
87 */
88
89 /* this routine used to be called 'rcomp' (WDP) */
90 /*---------------------------------------------------------------------------*/
91
fits_rcomp(int a[],int nx,unsigned char * c,int clen,int nblock)92 int fits_rcomp(int a[], /* input array */
93 int nx, /* number of input pixels */
94 unsigned char *c, /* output buffer */
95 int clen, /* max length of output */
96 int nblock) /* coding block size */
97 {
98 Buffer bufmem, *buffer = &bufmem;
99 /* int bsize; */
100 int i, j, thisblock;
101 int lastpix, nextpix, pdiff;
102 int v, fs, fsmask, top, fsmax, fsbits, bbits;
103 int lbitbuffer, lbits_to_go;
104 unsigned int psum;
105 double pixelsum, dpsum;
106 unsigned int *diff;
107
108 /*
109 * Original size of each pixel (bsize, bytes) and coding block
110 * size (nblock, pixels)
111 * Could make bsize a parameter to allow more efficient
112 * compression of short & byte images.
113 */
114 /* bsize = 4; */
115
116 /* nblock = 32; now an input parameter*/
117 /*
118 * From bsize derive:
119 * FSBITS = # bits required to store FS
120 * FSMAX = maximum value for FS
121 * BBITS = bits/pixel for direct coding
122 */
123
124 /*
125 switch (bsize) {
126 case 1:
127 fsbits = 3;
128 fsmax = 6;
129 break;
130 case 2:
131 fsbits = 4;
132 fsmax = 14;
133 break;
134 case 4:
135 fsbits = 5;
136 fsmax = 25;
137 break;
138 default:
139 ffpmsg("rdecomp: bsize must be 1, 2, or 4 bytes");
140 return(-1);
141 }
142 */
143
144 /* move out of switch block, to tweak performance */
145 fsbits = 5;
146 fsmax = 25;
147
148 bbits = 1<<fsbits;
149
150 /*
151 * Set up buffer pointers
152 */
153 buffer->start = c;
154 buffer->current = c;
155 buffer->end = c+clen;
156 buffer->bits_to_go = 8;
157 /*
158 * array for differences mapped to non-negative values
159 */
160 diff = (unsigned int *) malloc(nblock*sizeof(unsigned int));
161 if (diff == (unsigned int *) NULL) {
162 ffpmsg("fits_rcomp: insufficient memory");
163 return(-1);
164 }
165 /*
166 * Code in blocks of nblock pixels
167 */
168 start_outputing_bits(buffer);
169
170 /* write out first int value to the first 4 bytes of the buffer */
171 if (output_nbits(buffer, a[0], 32) == EOF) {
172 ffpmsg("rice_encode: end of buffer");
173 free(diff);
174 return(-1);
175 }
176
177 lastpix = a[0]; /* the first difference will always be zero */
178
179 thisblock = nblock;
180 for (i=0; i<nx; i += nblock) {
181 /* last block may be shorter */
182 if (nx-i < nblock) thisblock = nx-i;
183 /*
184 * Compute differences of adjacent pixels and map them to unsigned values.
185 * Note that this may overflow the integer variables -- that's
186 * OK, because we can recover when decompressing. If we were
187 * compressing shorts or bytes, would want to do this arithmetic
188 * with short/byte working variables (though diff will still be
189 * passed as an int.)
190 *
191 * compute sum of mapped pixel values at same time
192 * use double precision for sum to allow 32-bit integer inputs
193 */
194 pixelsum = 0.0;
195 for (j=0; j<thisblock; j++) {
196 nextpix = a[i+j];
197 pdiff = nextpix - lastpix;
198 diff[j] = (unsigned int) ((pdiff<0) ? ~(pdiff<<1) : (pdiff<<1));
199 pixelsum += diff[j];
200 lastpix = nextpix;
201 }
202
203 /*
204 * compute number of bits to split from sum
205 */
206 dpsum = (pixelsum - (thisblock/2) - 1)/thisblock;
207 if (dpsum < 0) dpsum = 0.0;
208 psum = ((unsigned int) dpsum ) >> 1;
209 for (fs = 0; psum>0; fs++) psum >>= 1;
210
211 /*
212 * write the codes
213 * fsbits ID bits used to indicate split level
214 */
215 if (fs >= fsmax) {
216 /* Special high entropy case when FS >= fsmax
217 * Just write pixel difference values directly, no Rice coding at all.
218 */
219 if (output_nbits(buffer, fsmax+1, fsbits) == EOF) {
220 ffpmsg("rice_encode: end of buffer");
221 free(diff);
222 return(-1);
223 }
224 for (j=0; j<thisblock; j++) {
225 if (output_nbits(buffer, diff[j], bbits) == EOF) {
226 ffpmsg("rice_encode: end of buffer");
227 free(diff);
228 return(-1);
229 }
230 }
231 } else if (fs == 0 && pixelsum == 0) {
232 /*
233 * special low entropy case when FS = 0 and pixelsum=0 (all
234 * pixels in block are zero.)
235 * Output a 0 and return
236 */
237 if (output_nbits(buffer, 0, fsbits) == EOF) {
238 ffpmsg("rice_encode: end of buffer");
239 free(diff);
240 return(-1);
241 }
242 } else {
243 /* normal case: not either very high or very low entropy */
244 if (output_nbits(buffer, fs+1, fsbits) == EOF) {
245 ffpmsg("rice_encode: end of buffer");
246 free(diff);
247 return(-1);
248 }
249 fsmask = (1<<fs) - 1;
250 /*
251 * local copies of bit buffer to improve optimization
252 */
253 lbitbuffer = buffer->bitbuffer;
254 lbits_to_go = buffer->bits_to_go;
255 for (j=0; j<thisblock; j++) {
256 v = diff[j];
257 top = v >> fs;
258 /*
259 * top is coded by top zeros + 1
260 */
261 if (lbits_to_go >= top+1) {
262 lbitbuffer <<= top+1;
263 lbitbuffer |= 1;
264 lbits_to_go -= top+1;
265 } else {
266 lbitbuffer <<= lbits_to_go;
267 putcbuf(lbitbuffer & 0xff,buffer);
268
269 for (top -= lbits_to_go; top>=8; top -= 8) {
270 putcbuf(0, buffer);
271 }
272 lbitbuffer = 1;
273 lbits_to_go = 7-top;
274 }
275 /*
276 * bottom FS bits are written without coding
277 * code is output_nbits, moved into this routine to reduce overheads
278 * This code potentially breaks if FS>24, so I am limiting
279 * FS to 24 by choice of FSMAX above.
280 */
281 if (fs > 0) {
282 lbitbuffer <<= fs;
283 lbitbuffer |= v & fsmask;
284 lbits_to_go -= fs;
285 while (lbits_to_go <= 0) {
286 putcbuf((lbitbuffer>>(-lbits_to_go)) & 0xff,buffer);
287 lbits_to_go += 8;
288 }
289 }
290 }
291
292 /* check if overflowed output buffer */
293 if (buffer->current > buffer->end) {
294 ffpmsg("rice_encode: end of buffer");
295 free(diff);
296 return(-1);
297 }
298 buffer->bitbuffer = lbitbuffer;
299 buffer->bits_to_go = lbits_to_go;
300 }
301 }
302 done_outputing_bits(buffer);
303 free(diff);
304 /*
305 * return number of bytes used
306 */
307 return(buffer->current - buffer->start);
308 }
309 /*---------------------------------------------------------------------------*/
310
fits_rcomp_short(short a[],int nx,unsigned char * c,int clen,int nblock)311 int fits_rcomp_short(
312 short a[], /* input array */
313 int nx, /* number of input pixels */
314 unsigned char *c, /* output buffer */
315 int clen, /* max length of output */
316 int nblock) /* coding block size */
317 {
318 Buffer bufmem, *buffer = &bufmem;
319 /* int bsize; */
320 int i, j, thisblock;
321
322 /*
323 NOTE: in principle, the following 2 variable could be declared as 'short'
324 but in fact the code runs faster (on 32-bit Linux at least) as 'int'
325 */
326 int lastpix, nextpix;
327 /* int pdiff; */
328 short pdiff;
329 int v, fs, fsmask, top, fsmax, fsbits, bbits;
330 int lbitbuffer, lbits_to_go;
331 /* unsigned int psum; */
332 unsigned short psum;
333 double pixelsum, dpsum;
334 unsigned int *diff;
335
336 /*
337 * Original size of each pixel (bsize, bytes) and coding block
338 * size (nblock, pixels)
339 * Could make bsize a parameter to allow more efficient
340 * compression of short & byte images.
341 */
342 /* bsize = 2; */
343
344 /* nblock = 32; now an input parameter */
345 /*
346 * From bsize derive:
347 * FSBITS = # bits required to store FS
348 * FSMAX = maximum value for FS
349 * BBITS = bits/pixel for direct coding
350 */
351
352 /*
353 switch (bsize) {
354 case 1:
355 fsbits = 3;
356 fsmax = 6;
357 break;
358 case 2:
359 fsbits = 4;
360 fsmax = 14;
361 break;
362 case 4:
363 fsbits = 5;
364 fsmax = 25;
365 break;
366 default:
367 ffpmsg("rdecomp: bsize must be 1, 2, or 4 bytes");
368 return(-1);
369 }
370 */
371
372 /* move these out of switch block to further tweak performance */
373 fsbits = 4;
374 fsmax = 14;
375
376 bbits = 1<<fsbits;
377
378 /*
379 * Set up buffer pointers
380 */
381 buffer->start = c;
382 buffer->current = c;
383 buffer->end = c+clen;
384 buffer->bits_to_go = 8;
385 /*
386 * array for differences mapped to non-negative values
387 */
388 diff = (unsigned int *) malloc(nblock*sizeof(unsigned int));
389 if (diff == (unsigned int *) NULL) {
390 ffpmsg("fits_rcomp: insufficient memory");
391 return(-1);
392 }
393 /*
394 * Code in blocks of nblock pixels
395 */
396 start_outputing_bits(buffer);
397
398 /* write out first short value to the first 2 bytes of the buffer */
399 if (output_nbits(buffer, a[0], 16) == EOF) {
400 ffpmsg("rice_encode: end of buffer");
401 free(diff);
402 return(-1);
403 }
404
405 lastpix = a[0]; /* the first difference will always be zero */
406
407 thisblock = nblock;
408 for (i=0; i<nx; i += nblock) {
409 /* last block may be shorter */
410 if (nx-i < nblock) thisblock = nx-i;
411 /*
412 * Compute differences of adjacent pixels and map them to unsigned values.
413 * Note that this may overflow the integer variables -- that's
414 * OK, because we can recover when decompressing. If we were
415 * compressing shorts or bytes, would want to do this arithmetic
416 * with short/byte working variables (though diff will still be
417 * passed as an int.)
418 *
419 * compute sum of mapped pixel values at same time
420 * use double precision for sum to allow 32-bit integer inputs
421 */
422 pixelsum = 0.0;
423 for (j=0; j<thisblock; j++) {
424 nextpix = a[i+j];
425 pdiff = nextpix - lastpix;
426 diff[j] = (unsigned int) ((pdiff<0) ? ~(pdiff<<1) : (pdiff<<1));
427 pixelsum += diff[j];
428 lastpix = nextpix;
429 }
430 /*
431 * compute number of bits to split from sum
432 */
433 dpsum = (pixelsum - (thisblock/2) - 1)/thisblock;
434 if (dpsum < 0) dpsum = 0.0;
435 /* psum = ((unsigned int) dpsum ) >> 1; */
436 psum = ((unsigned short) dpsum ) >> 1;
437 for (fs = 0; psum>0; fs++) psum >>= 1;
438
439 /*
440 * write the codes
441 * fsbits ID bits used to indicate split level
442 */
443 if (fs >= fsmax) {
444 /* case3++; */
445 /* Special high entropy case when FS >= fsmax
446 * Just write pixel difference values directly, no Rice coding at all.
447 */
448 if (output_nbits(buffer, fsmax+1, fsbits) == EOF) {
449 ffpmsg("rice_encode: end of buffer");
450 free(diff);
451 return(-1);
452 }
453 for (j=0; j<thisblock; j++) {
454 if (output_nbits(buffer, diff[j], bbits) == EOF) {
455 ffpmsg("rice_encode: end of buffer");
456 free(diff);
457 return(-1);
458 }
459 }
460 } else if (fs == 0 && pixelsum == 0) {
461 /* case1++; */
462 /*
463 * special low entropy case when FS = 0 and pixelsum=0 (all
464 * pixels in block are zero.)
465 * Output a 0 and return
466 */
467 if (output_nbits(buffer, 0, fsbits) == EOF) {
468 ffpmsg("rice_encode: end of buffer");
469 free(diff);
470 return(-1);
471 }
472 } else {
473 /* case2++; */
474 /* normal case: not either very high or very low entropy */
475 if (output_nbits(buffer, fs+1, fsbits) == EOF) {
476 ffpmsg("rice_encode: end of buffer");
477 free(diff);
478 return(-1);
479 }
480 fsmask = (1<<fs) - 1;
481 /*
482 * local copies of bit buffer to improve optimization
483 */
484 lbitbuffer = buffer->bitbuffer;
485 lbits_to_go = buffer->bits_to_go;
486 for (j=0; j<thisblock; j++) {
487 v = diff[j];
488 top = v >> fs;
489 /*
490 * top is coded by top zeros + 1
491 */
492 if (lbits_to_go >= top+1) {
493 lbitbuffer <<= top+1;
494 lbitbuffer |= 1;
495 lbits_to_go -= top+1;
496 } else {
497 lbitbuffer <<= lbits_to_go;
498 putcbuf(lbitbuffer & 0xff,buffer);
499 for (top -= lbits_to_go; top>=8; top -= 8) {
500 putcbuf(0, buffer);
501 }
502 lbitbuffer = 1;
503 lbits_to_go = 7-top;
504 }
505 /*
506 * bottom FS bits are written without coding
507 * code is output_nbits, moved into this routine to reduce overheads
508 * This code potentially breaks if FS>24, so I am limiting
509 * FS to 24 by choice of FSMAX above.
510 */
511 if (fs > 0) {
512 lbitbuffer <<= fs;
513 lbitbuffer |= v & fsmask;
514 lbits_to_go -= fs;
515 while (lbits_to_go <= 0) {
516 putcbuf((lbitbuffer>>(-lbits_to_go)) & 0xff,buffer);
517 lbits_to_go += 8;
518 }
519 }
520 }
521 /* check if overflowed output buffer */
522 if (buffer->current > buffer->end) {
523 ffpmsg("rice_encode: end of buffer");
524 free(diff);
525 return(-1);
526 }
527 buffer->bitbuffer = lbitbuffer;
528 buffer->bits_to_go = lbits_to_go;
529 }
530 }
531 done_outputing_bits(buffer);
532 free(diff);
533 /*
534 * return number of bytes used
535 */
536 return(buffer->current - buffer->start);
537 }
538 /*---------------------------------------------------------------------------*/
539
fits_rcomp_byte(signed char a[],int nx,unsigned char * c,int clen,int nblock)540 int fits_rcomp_byte(
541 signed char a[], /* input array */
542 int nx, /* number of input pixels */
543 unsigned char *c, /* output buffer */
544 int clen, /* max length of output */
545 int nblock) /* coding block size */
546 {
547 Buffer bufmem, *buffer = &bufmem;
548 /* int bsize; */
549 int i, j, thisblock;
550
551 /*
552 NOTE: in principle, the following 2 variable could be declared as 'short'
553 but in fact the code runs faster (on 32-bit Linux at least) as 'int'
554 */
555 int lastpix, nextpix;
556 /* int pdiff; */
557 signed char pdiff;
558 int v, fs, fsmask, top, fsmax, fsbits, bbits;
559 int lbitbuffer, lbits_to_go;
560 /* unsigned int psum; */
561 unsigned char psum;
562 double pixelsum, dpsum;
563 unsigned int *diff;
564
565 /*
566 * Original size of each pixel (bsize, bytes) and coding block
567 * size (nblock, pixels)
568 * Could make bsize a parameter to allow more efficient
569 * compression of short & byte images.
570 */
571 /* bsize = 1; */
572
573 /* nblock = 32; now an input parameter */
574 /*
575 * From bsize derive:
576 * FSBITS = # bits required to store FS
577 * FSMAX = maximum value for FS
578 * BBITS = bits/pixel for direct coding
579 */
580
581 /*
582 switch (bsize) {
583 case 1:
584 fsbits = 3;
585 fsmax = 6;
586 break;
587 case 2:
588 fsbits = 4;
589 fsmax = 14;
590 break;
591 case 4:
592 fsbits = 5;
593 fsmax = 25;
594 break;
595 default:
596 ffpmsg("rdecomp: bsize must be 1, 2, or 4 bytes");
597 return(-1);
598 }
599 */
600
601 /* move these out of switch block to further tweak performance */
602 fsbits = 3;
603 fsmax = 6;
604 bbits = 1<<fsbits;
605
606 /*
607 * Set up buffer pointers
608 */
609 buffer->start = c;
610 buffer->current = c;
611 buffer->end = c+clen;
612 buffer->bits_to_go = 8;
613 /*
614 * array for differences mapped to non-negative values
615 */
616 diff = (unsigned int *) malloc(nblock*sizeof(unsigned int));
617 if (diff == (unsigned int *) NULL) {
618 ffpmsg("fits_rcomp: insufficient memory");
619 return(-1);
620 }
621 /*
622 * Code in blocks of nblock pixels
623 */
624 start_outputing_bits(buffer);
625
626 /* write out first byte value to the first byte of the buffer */
627 if (output_nbits(buffer, a[0], 8) == EOF) {
628 ffpmsg("rice_encode: end of buffer");
629 free(diff);
630 return(-1);
631 }
632
633 lastpix = a[0]; /* the first difference will always be zero */
634
635 thisblock = nblock;
636 for (i=0; i<nx; i += nblock) {
637 /* last block may be shorter */
638 if (nx-i < nblock) thisblock = nx-i;
639 /*
640 * Compute differences of adjacent pixels and map them to unsigned values.
641 * Note that this may overflow the integer variables -- that's
642 * OK, because we can recover when decompressing. If we were
643 * compressing shorts or bytes, would want to do this arithmetic
644 * with short/byte working variables (though diff will still be
645 * passed as an int.)
646 *
647 * compute sum of mapped pixel values at same time
648 * use double precision for sum to allow 32-bit integer inputs
649 */
650 pixelsum = 0.0;
651 for (j=0; j<thisblock; j++) {
652 nextpix = a[i+j];
653 pdiff = nextpix - lastpix;
654 diff[j] = (unsigned int) ((pdiff<0) ? ~(pdiff<<1) : (pdiff<<1));
655 pixelsum += diff[j];
656 lastpix = nextpix;
657 }
658 /*
659 * compute number of bits to split from sum
660 */
661 dpsum = (pixelsum - (thisblock/2) - 1)/thisblock;
662 if (dpsum < 0) dpsum = 0.0;
663 /* psum = ((unsigned int) dpsum ) >> 1; */
664 psum = ((unsigned char) dpsum ) >> 1;
665 for (fs = 0; psum>0; fs++) psum >>= 1;
666
667 /*
668 * write the codes
669 * fsbits ID bits used to indicate split level
670 */
671 if (fs >= fsmax) {
672 /* Special high entropy case when FS >= fsmax
673 * Just write pixel difference values directly, no Rice coding at all.
674 */
675 if (output_nbits(buffer, fsmax+1, fsbits) == EOF) {
676 ffpmsg("rice_encode: end of buffer");
677 free(diff);
678 return(-1);
679 }
680 for (j=0; j<thisblock; j++) {
681 if (output_nbits(buffer, diff[j], bbits) == EOF) {
682 ffpmsg("rice_encode: end of buffer");
683 free(diff);
684 return(-1);
685 }
686 }
687 } else if (fs == 0 && pixelsum == 0) {
688 /*
689 * special low entropy case when FS = 0 and pixelsum=0 (all
690 * pixels in block are zero.)
691 * Output a 0 and return
692 */
693 if (output_nbits(buffer, 0, fsbits) == EOF) {
694 ffpmsg("rice_encode: end of buffer");
695 free(diff);
696 return(-1);
697 }
698 } else {
699 /* normal case: not either very high or very low entropy */
700 if (output_nbits(buffer, fs+1, fsbits) == EOF) {
701 ffpmsg("rice_encode: end of buffer");
702 free(diff);
703 return(-1);
704 }
705 fsmask = (1<<fs) - 1;
706 /*
707 * local copies of bit buffer to improve optimization
708 */
709 lbitbuffer = buffer->bitbuffer;
710 lbits_to_go = buffer->bits_to_go;
711 for (j=0; j<thisblock; j++) {
712 v = diff[j];
713 top = v >> fs;
714 /*
715 * top is coded by top zeros + 1
716 */
717 if (lbits_to_go >= top+1) {
718 lbitbuffer <<= top+1;
719 lbitbuffer |= 1;
720 lbits_to_go -= top+1;
721 } else {
722 lbitbuffer <<= lbits_to_go;
723 putcbuf(lbitbuffer & 0xff,buffer);
724 for (top -= lbits_to_go; top>=8; top -= 8) {
725 putcbuf(0, buffer);
726 }
727 lbitbuffer = 1;
728 lbits_to_go = 7-top;
729 }
730 /*
731 * bottom FS bits are written without coding
732 * code is output_nbits, moved into this routine to reduce overheads
733 * This code potentially breaks if FS>24, so I am limiting
734 * FS to 24 by choice of FSMAX above.
735 */
736 if (fs > 0) {
737 lbitbuffer <<= fs;
738 lbitbuffer |= v & fsmask;
739 lbits_to_go -= fs;
740 while (lbits_to_go <= 0) {
741 putcbuf((lbitbuffer>>(-lbits_to_go)) & 0xff,buffer);
742 lbits_to_go += 8;
743 }
744 }
745 }
746 /* check if overflowed output buffer */
747 if (buffer->current > buffer->end) {
748 ffpmsg("rice_encode: end of buffer");
749 free(diff);
750 return(-1);
751 }
752 buffer->bitbuffer = lbitbuffer;
753 buffer->bits_to_go = lbits_to_go;
754 }
755 }
756 done_outputing_bits(buffer);
757 free(diff);
758 /*
759 * return number of bytes used
760 */
761 return(buffer->current - buffer->start);
762 }
763 /*---------------------------------------------------------------------------*/
764 /* bit_output.c
765 *
766 * Bit output routines
767 * Procedures return zero on success, EOF on end-of-buffer
768 *
769 * Programmer: R. White Date: 20 July 1998
770 */
771
772 /* Initialize for bit output */
773
start_outputing_bits(Buffer * buffer)774 static void start_outputing_bits(Buffer *buffer)
775 {
776 /*
777 * Buffer is empty to start with
778 */
779 buffer->bitbuffer = 0;
780 buffer->bits_to_go = 8;
781 }
782
783 /*---------------------------------------------------------------------------*/
784 /* Output N bits (N must be <= 32) */
785
output_nbits(Buffer * buffer,int bits,int n)786 static int output_nbits(Buffer *buffer, int bits, int n)
787 {
788 /* local copies */
789 int lbitbuffer;
790 int lbits_to_go;
791 /* AND mask for the right-most n bits */
792 static unsigned int mask[33] =
793 {0,
794 0x1, 0x3, 0x7, 0xf, 0x1f, 0x3f, 0x7f, 0xff,
795 0x1ff, 0x3ff, 0x7ff, 0xfff, 0x1fff, 0x3fff, 0x7fff, 0xffff,
796 0x1ffff, 0x3ffff, 0x7ffff, 0xfffff, 0x1fffff, 0x3fffff, 0x7fffff, 0xffffff,
797 0x1ffffff, 0x3ffffff, 0x7ffffff, 0xfffffff, 0x1fffffff, 0x3fffffff, 0x7fffffff, 0xffffffff};
798
799 /*
800 * insert bits at end of bitbuffer
801 */
802 lbitbuffer = buffer->bitbuffer;
803 lbits_to_go = buffer->bits_to_go;
804 if (lbits_to_go+n > 32) {
805 /*
806 * special case for large n: put out the top lbits_to_go bits first
807 * note that 0 < lbits_to_go <= 8
808 */
809 lbitbuffer <<= lbits_to_go;
810 /* lbitbuffer |= (bits>>(n-lbits_to_go)) & ((1<<lbits_to_go)-1); */
811 lbitbuffer |= (bits>>(n-lbits_to_go)) & *(mask+lbits_to_go);
812 putcbuf(lbitbuffer & 0xff,buffer);
813 n -= lbits_to_go;
814 lbits_to_go = 8;
815 }
816 lbitbuffer <<= n;
817 /* lbitbuffer |= ( bits & ((1<<n)-1) ); */
818 lbitbuffer |= ( bits & *(mask+n) );
819 lbits_to_go -= n;
820 while (lbits_to_go <= 0) {
821 /*
822 * bitbuffer full, put out top 8 bits
823 */
824 putcbuf((lbitbuffer>>(-lbits_to_go)) & 0xff,buffer);
825 lbits_to_go += 8;
826 }
827 buffer->bitbuffer = lbitbuffer;
828 buffer->bits_to_go = lbits_to_go;
829 return(0);
830 }
831 /*---------------------------------------------------------------------------*/
832 /* Flush out the last bits */
833
done_outputing_bits(Buffer * buffer)834 static int done_outputing_bits(Buffer *buffer)
835 {
836 if(buffer->bits_to_go < 8) {
837 putcbuf(buffer->bitbuffer<<buffer->bits_to_go,buffer);
838
839 /* if (putcbuf(buffer->bitbuffer<<buffer->bits_to_go,buffer) == EOF)
840 return(EOF);
841 */
842 }
843 return(0);
844 }
845 /*---------------------------------------------------------------------------*/
846 /*----------------------------------------------------------*/
847 /* */
848 /* START OF SOURCE FILE ORIGINALLY CALLED rdecomp.c */
849 /* */
850 /*----------------------------------------------------------*/
851
852 /* @(#) rdecomp.c 1.4 99/03/01 12:38:41 */
853 /* rdecomp.c Decompress image line using
854 * (1) Difference of adjacent pixels
855 * (2) Rice algorithm coding
856 *
857 * Returns 0 on success or 1 on failure
858 */
859
860 /* moved these 'includes' to the beginning of the file (WDP)
861 #include <stdio.h>
862 #include <stdlib.h>
863 */
864
865 /*---------------------------------------------------------------------------*/
866 /* this routine used to be called 'rdecomp' (WDP) */
867
fits_rdecomp(unsigned char * c,int clen,unsigned int array[],int nx,int nblock)868 int fits_rdecomp (unsigned char *c, /* input buffer */
869 int clen, /* length of input */
870 unsigned int array[], /* output array */
871 int nx, /* number of output pixels */
872 int nblock) /* coding block size */
873 {
874 /* int bsize; */
875 int i, k, imax;
876 int nbits, nzero, fs;
877 unsigned char *cend, bytevalue;
878 unsigned int b, diff, lastpix;
879 int fsmax, fsbits, bbits;
880 extern const int nonzero_count[];
881
882 /*
883 * Original size of each pixel (bsize, bytes) and coding block
884 * size (nblock, pixels)
885 * Could make bsize a parameter to allow more efficient
886 * compression of short & byte images.
887 */
888 /* bsize = 4; */
889
890 /* nblock = 32; now an input parameter */
891 /*
892 * From bsize derive:
893 * FSBITS = # bits required to store FS
894 * FSMAX = maximum value for FS
895 * BBITS = bits/pixel for direct coding
896 */
897
898 /*
899 switch (bsize) {
900 case 1:
901 fsbits = 3;
902 fsmax = 6;
903 break;
904 case 2:
905 fsbits = 4;
906 fsmax = 14;
907 break;
908 case 4:
909 fsbits = 5;
910 fsmax = 25;
911 break;
912 default:
913 ffpmsg("rdecomp: bsize must be 1, 2, or 4 bytes");
914 return 1;
915 }
916 */
917
918 /* move out of switch block, to tweak performance */
919 fsbits = 5;
920 fsmax = 25;
921
922 bbits = 1<<fsbits;
923
924 /*
925 * Decode in blocks of nblock pixels
926 */
927
928 /* first 4 bytes of input buffer contain the value of the first */
929 /* 4 byte integer value, without any encoding */
930
931 lastpix = 0;
932 bytevalue = c[0];
933 lastpix = lastpix | (bytevalue<<24);
934 bytevalue = c[1];
935 lastpix = lastpix | (bytevalue<<16);
936 bytevalue = c[2];
937 lastpix = lastpix | (bytevalue<<8);
938 bytevalue = c[3];
939 lastpix = lastpix | bytevalue;
940
941 c += 4;
942 cend = c + clen - 4;
943
944 b = *c++; /* bit buffer */
945 nbits = 8; /* number of bits remaining in b */
946 for (i = 0; i<nx; ) {
947 /* get the FS value from first fsbits */
948 nbits -= fsbits;
949 while (nbits < 0) {
950 b = (b<<8) | (*c++);
951 nbits += 8;
952 }
953 fs = (b >> nbits) - 1;
954
955 b &= (1<<nbits)-1;
956 /* loop over the next block */
957 imax = i + nblock;
958 if (imax > nx) imax = nx;
959 if (fs<0) {
960 /* low-entropy case, all zero differences */
961 for ( ; i<imax; i++) array[i] = lastpix;
962 } else if (fs==fsmax) {
963 /* high-entropy case, directly coded pixel values */
964 for ( ; i<imax; i++) {
965 k = bbits - nbits;
966 diff = b<<k;
967 for (k -= 8; k >= 0; k -= 8) {
968 b = *c++;
969 diff |= b<<k;
970 }
971 if (nbits>0) {
972 b = *c++;
973 diff |= b>>(-k);
974 b &= (1<<nbits)-1;
975 } else {
976 b = 0;
977 }
978 /*
979 * undo mapping and differencing
980 * Note that some of these operations will overflow the
981 * unsigned int arithmetic -- that's OK, it all works
982 * out to give the right answers in the output file.
983 */
984 if ((diff & 1) == 0) {
985 diff = diff>>1;
986 } else {
987 diff = ~(diff>>1);
988 }
989 array[i] = diff+lastpix;
990 lastpix = array[i];
991 }
992 } else {
993 /* normal case, Rice coding */
994 for ( ; i<imax; i++) {
995 /* count number of leading zeros */
996 while (b == 0) {
997 nbits += 8;
998 b = *c++;
999 }
1000 nzero = nbits - nonzero_count[b];
1001 nbits -= nzero+1;
1002 /* flip the leading one-bit */
1003 b ^= 1<<nbits;
1004 /* get the FS trailing bits */
1005 nbits -= fs;
1006 while (nbits < 0) {
1007 b = (b<<8) | (*c++);
1008 nbits += 8;
1009 }
1010 diff = (nzero<<fs) | (b>>nbits);
1011 b &= (1<<nbits)-1;
1012
1013 /* undo mapping and differencing */
1014 if ((diff & 1) == 0) {
1015 diff = diff>>1;
1016 } else {
1017 diff = ~(diff>>1);
1018 }
1019 array[i] = diff+lastpix;
1020 lastpix = array[i];
1021 }
1022 }
1023 if (c > cend) {
1024 ffpmsg("decompression error: hit end of compressed byte stream");
1025 return 1;
1026 }
1027 }
1028 if (c < cend) {
1029 ffpmsg("decompression warning: unused bytes at end of compressed buffer");
1030 }
1031 return 0;
1032 }
1033 /*---------------------------------------------------------------------------*/
1034 /* this routine used to be called 'rdecomp' (WDP) */
1035
fits_rdecomp_short(unsigned char * c,int clen,unsigned short array[],int nx,int nblock)1036 int fits_rdecomp_short (unsigned char *c, /* input buffer */
1037 int clen, /* length of input */
1038 unsigned short array[], /* output array */
1039 int nx, /* number of output pixels */
1040 int nblock) /* coding block size */
1041 {
1042 int i, imax;
1043 /* int bsize; */
1044 int k;
1045 int nbits, nzero, fs;
1046 unsigned char *cend, bytevalue;
1047 unsigned int b, diff, lastpix;
1048 int fsmax, fsbits, bbits;
1049 extern const int nonzero_count[];
1050
1051 /*
1052 * Original size of each pixel (bsize, bytes) and coding block
1053 * size (nblock, pixels)
1054 * Could make bsize a parameter to allow more efficient
1055 * compression of short & byte images.
1056 */
1057
1058 /* bsize = 2; */
1059
1060 /* nblock = 32; now an input parameter */
1061 /*
1062 * From bsize derive:
1063 * FSBITS = # bits required to store FS
1064 * FSMAX = maximum value for FS
1065 * BBITS = bits/pixel for direct coding
1066 */
1067
1068 /*
1069 switch (bsize) {
1070 case 1:
1071 fsbits = 3;
1072 fsmax = 6;
1073 break;
1074 case 2:
1075 fsbits = 4;
1076 fsmax = 14;
1077 break;
1078 case 4:
1079 fsbits = 5;
1080 fsmax = 25;
1081 break;
1082 default:
1083 ffpmsg("rdecomp: bsize must be 1, 2, or 4 bytes");
1084 return 1;
1085 }
1086 */
1087
1088 /* move out of switch block, to tweak performance */
1089 fsbits = 4;
1090 fsmax = 14;
1091
1092 bbits = 1<<fsbits;
1093
1094 /*
1095 * Decode in blocks of nblock pixels
1096 */
1097
1098 /* first 2 bytes of input buffer contain the value of the first */
1099 /* 2 byte integer value, without any encoding */
1100
1101 lastpix = 0;
1102 bytevalue = c[0];
1103 lastpix = lastpix | (bytevalue<<8);
1104 bytevalue = c[1];
1105 lastpix = lastpix | bytevalue;
1106
1107 c += 2;
1108 cend = c + clen - 2;
1109
1110 b = *c++; /* bit buffer */
1111 nbits = 8; /* number of bits remaining in b */
1112 for (i = 0; i<nx; ) {
1113 /* get the FS value from first fsbits */
1114 nbits -= fsbits;
1115 while (nbits < 0) {
1116 b = (b<<8) | (*c++);
1117 nbits += 8;
1118 }
1119 fs = (b >> nbits) - 1;
1120
1121 b &= (1<<nbits)-1;
1122 /* loop over the next block */
1123 imax = i + nblock;
1124 if (imax > nx) imax = nx;
1125 if (fs<0) {
1126 /* low-entropy case, all zero differences */
1127 for ( ; i<imax; i++) array[i] = lastpix;
1128 } else if (fs==fsmax) {
1129 /* high-entropy case, directly coded pixel values */
1130 for ( ; i<imax; i++) {
1131 k = bbits - nbits;
1132 diff = b<<k;
1133 for (k -= 8; k >= 0; k -= 8) {
1134 b = *c++;
1135 diff |= b<<k;
1136 }
1137 if (nbits>0) {
1138 b = *c++;
1139 diff |= b>>(-k);
1140 b &= (1<<nbits)-1;
1141 } else {
1142 b = 0;
1143 }
1144
1145 /*
1146 * undo mapping and differencing
1147 * Note that some of these operations will overflow the
1148 * unsigned int arithmetic -- that's OK, it all works
1149 * out to give the right answers in the output file.
1150 */
1151 if ((diff & 1) == 0) {
1152 diff = diff>>1;
1153 } else {
1154 diff = ~(diff>>1);
1155 }
1156 array[i] = diff+lastpix;
1157 lastpix = array[i];
1158 }
1159 } else {
1160 /* normal case, Rice coding */
1161 for ( ; i<imax; i++) {
1162 /* count number of leading zeros */
1163 while (b == 0) {
1164 nbits += 8;
1165 b = *c++;
1166 }
1167 nzero = nbits - nonzero_count[b];
1168 nbits -= nzero+1;
1169 /* flip the leading one-bit */
1170 b ^= 1<<nbits;
1171 /* get the FS trailing bits */
1172 nbits -= fs;
1173 while (nbits < 0) {
1174 b = (b<<8) | (*c++);
1175 nbits += 8;
1176 }
1177 diff = (nzero<<fs) | (b>>nbits);
1178 b &= (1<<nbits)-1;
1179
1180 /* undo mapping and differencing */
1181 if ((diff & 1) == 0) {
1182 diff = diff>>1;
1183 } else {
1184 diff = ~(diff>>1);
1185 }
1186 array[i] = diff+lastpix;
1187 lastpix = array[i];
1188 }
1189 }
1190 if (c > cend) {
1191 ffpmsg("decompression error: hit end of compressed byte stream");
1192 return 1;
1193 }
1194 }
1195 if (c < cend) {
1196 ffpmsg("decompression warning: unused bytes at end of compressed buffer");
1197 }
1198 return 0;
1199 }
1200 /*---------------------------------------------------------------------------*/
1201 /* this routine used to be called 'rdecomp' (WDP) */
1202
fits_rdecomp_byte(unsigned char * c,int clen,unsigned char array[],int nx,int nblock)1203 int fits_rdecomp_byte (unsigned char *c, /* input buffer */
1204 int clen, /* length of input */
1205 unsigned char array[], /* output array */
1206 int nx, /* number of output pixels */
1207 int nblock) /* coding block size */
1208 {
1209 int i, imax;
1210 /* int bsize; */
1211 int k;
1212 int nbits, nzero, fs;
1213 unsigned char *cend;
1214 unsigned int b, diff, lastpix;
1215 int fsmax, fsbits, bbits;
1216 extern const int nonzero_count[];
1217
1218 /*
1219 * Original size of each pixel (bsize, bytes) and coding block
1220 * size (nblock, pixels)
1221 * Could make bsize a parameter to allow more efficient
1222 * compression of short & byte images.
1223 */
1224
1225 /* bsize = 1; */
1226
1227 /* nblock = 32; now an input parameter */
1228 /*
1229 * From bsize derive:
1230 * FSBITS = # bits required to store FS
1231 * FSMAX = maximum value for FS
1232 * BBITS = bits/pixel for direct coding
1233 */
1234
1235 /*
1236 switch (bsize) {
1237 case 1:
1238 fsbits = 3;
1239 fsmax = 6;
1240 break;
1241 case 2:
1242 fsbits = 4;
1243 fsmax = 14;
1244 break;
1245 case 4:
1246 fsbits = 5;
1247 fsmax = 25;
1248 break;
1249 default:
1250 ffpmsg("rdecomp: bsize must be 1, 2, or 4 bytes");
1251 return 1;
1252 }
1253 */
1254
1255 /* move out of switch block, to tweak performance */
1256 fsbits = 3;
1257 fsmax = 6;
1258
1259 bbits = 1<<fsbits;
1260
1261 /*
1262 * Decode in blocks of nblock pixels
1263 */
1264
1265 /* first byte of input buffer contain the value of the first */
1266 /* byte integer value, without any encoding */
1267
1268 lastpix = c[0];
1269 c += 1;
1270 cend = c + clen - 1;
1271
1272 b = *c++; /* bit buffer */
1273 nbits = 8; /* number of bits remaining in b */
1274 for (i = 0; i<nx; ) {
1275 /* get the FS value from first fsbits */
1276 nbits -= fsbits;
1277 while (nbits < 0) {
1278 b = (b<<8) | (*c++);
1279 nbits += 8;
1280 }
1281 fs = (b >> nbits) - 1;
1282
1283 b &= (1<<nbits)-1;
1284 /* loop over the next block */
1285 imax = i + nblock;
1286 if (imax > nx) imax = nx;
1287 if (fs<0) {
1288 /* low-entropy case, all zero differences */
1289 for ( ; i<imax; i++) array[i] = lastpix;
1290 } else if (fs==fsmax) {
1291 /* high-entropy case, directly coded pixel values */
1292 for ( ; i<imax; i++) {
1293 k = bbits - nbits;
1294 diff = b<<k;
1295 for (k -= 8; k >= 0; k -= 8) {
1296 b = *c++;
1297 diff |= b<<k;
1298 }
1299 if (nbits>0) {
1300 b = *c++;
1301 diff |= b>>(-k);
1302 b &= (1<<nbits)-1;
1303 } else {
1304 b = 0;
1305 }
1306
1307 /*
1308 * undo mapping and differencing
1309 * Note that some of these operations will overflow the
1310 * unsigned int arithmetic -- that's OK, it all works
1311 * out to give the right answers in the output file.
1312 */
1313 if ((diff & 1) == 0) {
1314 diff = diff>>1;
1315 } else {
1316 diff = ~(diff>>1);
1317 }
1318 array[i] = diff+lastpix;
1319 lastpix = array[i];
1320 }
1321 } else {
1322 /* normal case, Rice coding */
1323 for ( ; i<imax; i++) {
1324 /* count number of leading zeros */
1325 while (b == 0) {
1326 nbits += 8;
1327 b = *c++;
1328 }
1329 nzero = nbits - nonzero_count[b];
1330 nbits -= nzero+1;
1331 /* flip the leading one-bit */
1332 b ^= 1<<nbits;
1333 /* get the FS trailing bits */
1334 nbits -= fs;
1335 while (nbits < 0) {
1336 b = (b<<8) | (*c++);
1337 nbits += 8;
1338 }
1339 diff = (nzero<<fs) | (b>>nbits);
1340 b &= (1<<nbits)-1;
1341
1342 /* undo mapping and differencing */
1343 if ((diff & 1) == 0) {
1344 diff = diff>>1;
1345 } else {
1346 diff = ~(diff>>1);
1347 }
1348 array[i] = diff+lastpix;
1349 lastpix = array[i];
1350 }
1351 }
1352 if (c > cend) {
1353 ffpmsg("decompression error: hit end of compressed byte stream");
1354 return 1;
1355 }
1356 }
1357 if (c < cend) {
1358 ffpmsg("decompression warning: unused bytes at end of compressed buffer");
1359 }
1360 return 0;
1361 }
1362