- Timestamp:
- Jul 3, 2008, 10:26:51 PM (18 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/eam_branch_20080421/Ohana/src/libfits/extern/ricecomp.c
r15487 r18424 1 2 // 2008.06.22 EAM : This code was taken from CFITSIO and included in Ohana for rice 3 // decompression. In order to include this .c file (or future upgrades), we include our own 4 // version of ricecomp.h (not the one available in the CFITSIO tree) and include it instead of 5 // fitsio2.h. This file defines our own version of ffpmsg. 6 1 7 /* 2 8 The following code was written by Richard White at STScI and made … … 23 29 #include <stdlib.h> 24 30 #include <string.h> 25 #include "ricecomp.h" /* originally included in rcomp.c file (WDP) */ 31 32 typedef unsigned char Buffer_t; 33 34 typedef struct { 35 int bitbuffer; /* bit buffer */ 36 int bits_to_go; /* bits to go in buffer */ 37 Buffer_t *start; /* start of buffer */ 38 Buffer_t *current; /* current position in buffer */ 39 Buffer_t *end; /* end of buffer */ 40 } Buffer; 41 42 #define putcbuf(c,mf) ((*(mf->current)++ = c), 0) 43 26 44 // #include "fitsio2.h" 27 45 # include "ricecomp.h" 28 46 29 47 static void start_outputing_bits(Buffer *buffer); … … 32 50 33 51 /* this routine used to be called 'rcomp' (WDP) */ 52 /*---------------------------------------------------------------------------*/ 34 53 35 54 int fits_rcomp(int a[], /* input array */ … … 55 74 */ 56 75 bsize = 4; 57 /* nblock = 32; */ 76 77 /* nblock = 32; now an input parameter*/ 58 78 /* 59 79 * From bsize derive: … … 62 82 * BBITS = bits/pixel for direct coding 63 83 */ 84 85 /* 64 86 switch (bsize) { 65 87 case 1: … … 76 98 break; 77 99 default: 78 f printf (stderr,"rdecomp: bsize must be 1, 2, or 4 bytes");100 ffpmsg("rdecomp: bsize must be 1, 2, or 4 bytes"); 79 101 return(-1); 80 102 } 103 */ 104 105 /* move out of switch block, to tweak performance */ 106 fsbits = 5; 107 fsmax = 25; 81 108 bbits = 1<<fsbits; 82 109 … … 93 120 diff = (unsigned int *) malloc(nblock*sizeof(unsigned int)); 94 121 if (diff == (unsigned int *) NULL) { 95 f printf (stderr,"fits_rcomp: insufficient memory");122 ffpmsg("fits_rcomp: insufficient memory"); 96 123 return(-1); 97 124 } … … 103 130 /* write out first int value to the first 4 bytes of the buffer */ 104 131 if (output_nbits(buffer, a[0], 32) == EOF) { 105 f printf (stderr,"rice_encode: end of buffer");132 ffpmsg("rice_encode: end of buffer"); 106 133 free(diff); 107 134 return(-1); … … 133 160 lastpix = nextpix; 134 161 } 162 135 163 /* 136 164 * compute number of bits to split from sum … … 140 168 psum = ((unsigned int) dpsum ) >> 1; 141 169 for (fs = 0; psum>0; fs++) psum >>= 1; 170 142 171 /* 143 172 * write the codes … … 149 178 */ 150 179 if (output_nbits(buffer, fsmax+1, fsbits) == EOF) { 151 f printf (stderr,"rice_encode: end of buffer");180 ffpmsg("rice_encode: end of buffer"); 152 181 free(diff); 153 182 return(-1); … … 155 184 for (j=0; j<thisblock; j++) { 156 185 if (output_nbits(buffer, diff[j], bbits) == EOF) { 157 f printf (stderr,"rice_encode: end of buffer");186 ffpmsg("rice_encode: end of buffer"); 158 187 free(diff); 159 188 return(-1); … … 167 196 */ 168 197 if (output_nbits(buffer, 0, fsbits) == EOF) { 169 f printf (stderr,"rice_encode: end of buffer");198 ffpmsg("rice_encode: end of buffer"); 170 199 free(diff); 171 200 return(-1); … … 174 203 /* normal case: not either very high or very low entropy */ 175 204 if (output_nbits(buffer, fs+1, fsbits) == EOF) { 176 f printf (stderr,"rice_encode: end of buffer");205 ffpmsg("rice_encode: end of buffer"); 177 206 free(diff); 178 207 return(-1); … … 196 225 } else { 197 226 lbitbuffer <<= lbits_to_go; 198 if (putcbuf(lbitbuffer & 0xff,buffer) == EOF) { 199 fprintf (stderr, "rice_encode: end of buffer"); 200 free(diff); 201 return(-1); 202 } 227 putcbuf(lbitbuffer & 0xff,buffer); 228 203 229 for (top -= lbits_to_go; top>=8; top -= 8) { 204 if (putcbuf(0, buffer) == EOF) { 205 fprintf (stderr, "rice_encode: end of buffer"); 206 free(diff); 207 return(-1); 208 } 230 putcbuf(0, buffer); 209 231 } 210 232 lbitbuffer = 1; … … 222 244 lbits_to_go -= fs; 223 245 while (lbits_to_go <= 0) { 224 if (putcbuf((lbitbuffer>>(-lbits_to_go)) & 0xff,buffer)==EOF) { 225 fprintf (stderr, "rice_encode: end of buffer"); 226 free(diff); 227 return(-1); 228 } 246 putcbuf((lbitbuffer>>(-lbits_to_go)) & 0xff,buffer); 229 247 lbits_to_go += 8; 230 248 } 231 249 } 250 } 251 252 /* check if overflowed output buffer */ 253 if (buffer->current > buffer->end) { 254 ffpmsg("rice_encode: end of buffer"); 255 free(diff); 256 return(-1); 257 } 258 buffer->bitbuffer = lbitbuffer; 259 buffer->bits_to_go = lbits_to_go; 260 } 261 } 262 done_outputing_bits(buffer); 263 free(diff); 264 /* 265 * return number of bytes used 266 */ 267 return(buffer->current - buffer->start); 268 } 269 /*---------------------------------------------------------------------------*/ 270 271 int fits_rcomp_short( 272 short a[], /* input array */ 273 int nx, /* number of input pixels */ 274 unsigned char *c, /* output buffer */ 275 int clen, /* max length of output */ 276 int nblock) /* coding block size */ 277 { 278 Buffer bufmem, *buffer = &bufmem; 279 int bsize, i, j, thisblock; 280 281 /* 282 NOTE: in principle, the following 2 variable could be declared as 'short' 283 but in fact the code runs faster (on 32-bit Linux at least) as 'int' 284 */ 285 int lastpix, nextpix; 286 /* int pdiff; */ 287 short pdiff; 288 int v, fs, fsmask, top, fsmax, fsbits, bbits; 289 int lbitbuffer, lbits_to_go; 290 /* unsigned int psum; */ 291 unsigned short psum; 292 double pixelsum, dpsum; 293 unsigned int *diff; 294 295 /* 296 * Original size of each pixel (bsize, bytes) and coding block 297 * size (nblock, pixels) 298 * Could make bsize a parameter to allow more efficient 299 * compression of short & byte images. 300 */ 301 bsize = 2; 302 303 /* nblock = 32; now an input parameter */ 304 /* 305 * From bsize derive: 306 * FSBITS = # bits required to store FS 307 * FSMAX = maximum value for FS 308 * BBITS = bits/pixel for direct coding 309 */ 310 311 /* 312 switch (bsize) { 313 case 1: 314 fsbits = 3; 315 fsmax = 6; 316 break; 317 case 2: 318 fsbits = 4; 319 fsmax = 14; 320 break; 321 case 4: 322 fsbits = 5; 323 fsmax = 25; 324 break; 325 default: 326 ffpmsg("rdecomp: bsize must be 1, 2, or 4 bytes"); 327 return(-1); 328 } 329 */ 330 331 /* move these out of switch block to further tweak performance */ 332 fsbits = 4; 333 fsmax = 14; 334 bbits = 1<<fsbits; 335 336 /* 337 * Set up buffer pointers 338 */ 339 buffer->start = c; 340 buffer->current = c; 341 buffer->end = c+clen; 342 buffer->bits_to_go = 8; 343 /* 344 * array for differences mapped to non-negative values 345 */ 346 diff = (unsigned int *) malloc(nblock*sizeof(unsigned int)); 347 if (diff == (unsigned int *) NULL) { 348 ffpmsg("fits_rcomp: insufficient memory"); 349 return(-1); 350 } 351 /* 352 * Code in blocks of nblock pixels 353 */ 354 start_outputing_bits(buffer); 355 356 /* write out first short value to the first 2 bytes of the buffer */ 357 if (output_nbits(buffer, a[0], 16) == EOF) { 358 ffpmsg("rice_encode: end of buffer"); 359 free(diff); 360 return(-1); 361 } 362 363 lastpix = a[0]; /* the first difference will always be zero */ 364 365 thisblock = nblock; 366 for (i=0; i<nx; i += nblock) { 367 /* last block may be shorter */ 368 if (nx-i < nblock) thisblock = nx-i; 369 /* 370 * Compute differences of adjacent pixels and map them to unsigned values. 371 * Note that this may overflow the integer variables -- that's 372 * OK, because we can recover when decompressing. If we were 373 * compressing shorts or bytes, would want to do this arithmetic 374 * with short/byte working variables (though diff will still be 375 * passed as an int.) 376 * 377 * compute sum of mapped pixel values at same time 378 * use double precision for sum to allow 32-bit integer inputs 379 */ 380 pixelsum = 0.0; 381 for (j=0; j<thisblock; j++) { 382 nextpix = a[i+j]; 383 pdiff = nextpix - lastpix; 384 diff[j] = (unsigned int) ((pdiff<0) ? ~(pdiff<<1) : (pdiff<<1)); 385 pixelsum += diff[j]; 386 lastpix = nextpix; 387 } 388 /* 389 * compute number of bits to split from sum 390 */ 391 dpsum = (pixelsum - (thisblock/2) - 1)/thisblock; 392 if (dpsum < 0) dpsum = 0.0; 393 /* psum = ((unsigned int) dpsum ) >> 1; */ 394 psum = ((unsigned short) dpsum ) >> 1; 395 for (fs = 0; psum>0; fs++) psum >>= 1; 396 397 /* 398 * write the codes 399 * fsbits ID bits used to indicate split level 400 */ 401 if (fs >= fsmax) { 402 /* Special high entropy case when FS >= fsmax 403 * Just write pixel difference values directly, no Rice coding at all. 404 */ 405 if (output_nbits(buffer, fsmax+1, fsbits) == EOF) { 406 ffpmsg("rice_encode: end of buffer"); 407 free(diff); 408 return(-1); 409 } 410 for (j=0; j<thisblock; j++) { 411 if (output_nbits(buffer, diff[j], bbits) == EOF) { 412 ffpmsg("rice_encode: end of buffer"); 413 free(diff); 414 return(-1); 415 } 416 } 417 } else if (fs == 0 && pixelsum == 0) { 418 /* 419 * special low entropy case when FS = 0 and pixelsum=0 (all 420 * pixels in block are zero.) 421 * Output a 0 and return 422 */ 423 if (output_nbits(buffer, 0, fsbits) == EOF) { 424 ffpmsg("rice_encode: end of buffer"); 425 free(diff); 426 return(-1); 427 } 428 } else { 429 /* normal case: not either very high or very low entropy */ 430 if (output_nbits(buffer, fs+1, fsbits) == EOF) { 431 ffpmsg("rice_encode: end of buffer"); 432 free(diff); 433 return(-1); 434 } 435 fsmask = (1<<fs) - 1; 436 /* 437 * local copies of bit buffer to improve optimization 438 */ 439 lbitbuffer = buffer->bitbuffer; 440 lbits_to_go = buffer->bits_to_go; 441 for (j=0; j<thisblock; j++) { 442 v = diff[j]; 443 top = v >> fs; 444 /* 445 * top is coded by top zeros + 1 446 */ 447 if (lbits_to_go >= top+1) { 448 lbitbuffer <<= top+1; 449 lbitbuffer |= 1; 450 lbits_to_go -= top+1; 451 } else { 452 lbitbuffer <<= lbits_to_go; 453 putcbuf(lbitbuffer & 0xff,buffer); 454 for (top -= lbits_to_go; top>=8; top -= 8) { 455 putcbuf(0, buffer); 456 } 457 lbitbuffer = 1; 458 lbits_to_go = 7-top; 459 } 460 /* 461 * bottom FS bits are written without coding 462 * code is output_nbits, moved into this routine to reduce overheads 463 * This code potentially breaks if FS>24, so I am limiting 464 * FS to 24 by choice of FSMAX above. 465 */ 466 if (fs > 0) { 467 lbitbuffer <<= fs; 468 lbitbuffer |= v & fsmask; 469 lbits_to_go -= fs; 470 while (lbits_to_go <= 0) { 471 putcbuf((lbitbuffer>>(-lbits_to_go)) & 0xff,buffer); 472 lbits_to_go += 8; 473 } 474 } 475 } 476 /* check if overflowed output buffer */ 477 if (buffer->current > buffer->end) { 478 ffpmsg("rice_encode: end of buffer"); 479 free(diff); 480 return(-1); 481 } 482 buffer->bitbuffer = lbitbuffer; 483 buffer->bits_to_go = lbits_to_go; 484 } 485 } 486 done_outputing_bits(buffer); 487 free(diff); 488 /* 489 * return number of bytes used 490 */ 491 return(buffer->current - buffer->start); 492 } 493 /*---------------------------------------------------------------------------*/ 494 495 int fits_rcomp_byte( 496 signed char a[], /* input array */ 497 int nx, /* number of input pixels */ 498 unsigned char *c, /* output buffer */ 499 int clen, /* max length of output */ 500 int nblock) /* coding block size */ 501 { 502 Buffer bufmem, *buffer = &bufmem; 503 int bsize, i, j, thisblock; 504 505 /* 506 NOTE: in principle, the following 2 variable could be declared as 'short' 507 but in fact the code runs faster (on 32-bit Linux at least) as 'int' 508 */ 509 int lastpix, nextpix; 510 /* int pdiff; */ 511 signed char pdiff; 512 int v, fs, fsmask, top, fsmax, fsbits, bbits; 513 int lbitbuffer, lbits_to_go; 514 /* unsigned int psum; */ 515 unsigned char psum; 516 double pixelsum, dpsum; 517 unsigned int *diff; 518 519 /* 520 * Original size of each pixel (bsize, bytes) and coding block 521 * size (nblock, pixels) 522 * Could make bsize a parameter to allow more efficient 523 * compression of short & byte images. 524 */ 525 bsize = 1; 526 527 /* nblock = 32; now an input parameter */ 528 /* 529 * From bsize derive: 530 * FSBITS = # bits required to store FS 531 * FSMAX = maximum value for FS 532 * BBITS = bits/pixel for direct coding 533 */ 534 535 /* 536 switch (bsize) { 537 case 1: 538 fsbits = 3; 539 fsmax = 6; 540 break; 541 case 2: 542 fsbits = 4; 543 fsmax = 14; 544 break; 545 case 4: 546 fsbits = 5; 547 fsmax = 25; 548 break; 549 default: 550 ffpmsg("rdecomp: bsize must be 1, 2, or 4 bytes"); 551 return(-1); 552 } 553 */ 554 555 /* move these out of switch block to further tweak performance */ 556 fsbits = 3; 557 fsmax = 6; 558 bbits = 1<<fsbits; 559 560 /* 561 * Set up buffer pointers 562 */ 563 buffer->start = c; 564 buffer->current = c; 565 buffer->end = c+clen; 566 buffer->bits_to_go = 8; 567 /* 568 * array for differences mapped to non-negative values 569 */ 570 diff = (unsigned int *) malloc(nblock*sizeof(unsigned int)); 571 if (diff == (unsigned int *) NULL) { 572 ffpmsg("fits_rcomp: insufficient memory"); 573 return(-1); 574 } 575 /* 576 * Code in blocks of nblock pixels 577 */ 578 start_outputing_bits(buffer); 579 580 /* write out first byte value to the first byte of the buffer */ 581 if (output_nbits(buffer, a[0], 8) == EOF) { 582 ffpmsg("rice_encode: end of buffer"); 583 free(diff); 584 return(-1); 585 } 586 587 lastpix = a[0]; /* the first difference will always be zero */ 588 589 thisblock = nblock; 590 for (i=0; i<nx; i += nblock) { 591 /* last block may be shorter */ 592 if (nx-i < nblock) thisblock = nx-i; 593 /* 594 * Compute differences of adjacent pixels and map them to unsigned values. 595 * Note that this may overflow the integer variables -- that's 596 * OK, because we can recover when decompressing. If we were 597 * compressing shorts or bytes, would want to do this arithmetic 598 * with short/byte working variables (though diff will still be 599 * passed as an int.) 600 * 601 * compute sum of mapped pixel values at same time 602 * use double precision for sum to allow 32-bit integer inputs 603 */ 604 pixelsum = 0.0; 605 for (j=0; j<thisblock; j++) { 606 nextpix = a[i+j]; 607 pdiff = nextpix - lastpix; 608 diff[j] = (unsigned int) ((pdiff<0) ? ~(pdiff<<1) : (pdiff<<1)); 609 pixelsum += diff[j]; 610 lastpix = nextpix; 611 } 612 /* 613 * compute number of bits to split from sum 614 */ 615 dpsum = (pixelsum - (thisblock/2) - 1)/thisblock; 616 if (dpsum < 0) dpsum = 0.0; 617 /* psum = ((unsigned int) dpsum ) >> 1; */ 618 psum = ((unsigned char) dpsum ) >> 1; 619 for (fs = 0; psum>0; fs++) psum >>= 1; 620 621 /* 622 * write the codes 623 * fsbits ID bits used to indicate split level 624 */ 625 if (fs >= fsmax) { 626 /* Special high entropy case when FS >= fsmax 627 * Just write pixel difference values directly, no Rice coding at all. 628 */ 629 if (output_nbits(buffer, fsmax+1, fsbits) == EOF) { 630 ffpmsg("rice_encode: end of buffer"); 631 free(diff); 632 return(-1); 633 } 634 for (j=0; j<thisblock; j++) { 635 if (output_nbits(buffer, diff[j], bbits) == EOF) { 636 ffpmsg("rice_encode: end of buffer"); 637 free(diff); 638 return(-1); 639 } 640 } 641 } else if (fs == 0 && pixelsum == 0) { 642 /* 643 * special low entropy case when FS = 0 and pixelsum=0 (all 644 * pixels in block are zero.) 645 * Output a 0 and return 646 */ 647 if (output_nbits(buffer, 0, fsbits) == EOF) { 648 ffpmsg("rice_encode: end of buffer"); 649 free(diff); 650 return(-1); 651 } 652 } else { 653 /* normal case: not either very high or very low entropy */ 654 if (output_nbits(buffer, fs+1, fsbits) == EOF) { 655 ffpmsg("rice_encode: end of buffer"); 656 free(diff); 657 return(-1); 658 } 659 fsmask = (1<<fs) - 1; 660 /* 661 * local copies of bit buffer to improve optimization 662 */ 663 lbitbuffer = buffer->bitbuffer; 664 lbits_to_go = buffer->bits_to_go; 665 for (j=0; j<thisblock; j++) { 666 v = diff[j]; 667 top = v >> fs; 668 /* 669 * top is coded by top zeros + 1 670 */ 671 if (lbits_to_go >= top+1) { 672 lbitbuffer <<= top+1; 673 lbitbuffer |= 1; 674 lbits_to_go -= top+1; 675 } else { 676 lbitbuffer <<= lbits_to_go; 677 putcbuf(lbitbuffer & 0xff,buffer); 678 for (top -= lbits_to_go; top>=8; top -= 8) { 679 putcbuf(0, buffer); 680 } 681 lbitbuffer = 1; 682 lbits_to_go = 7-top; 683 } 684 /* 685 * bottom FS bits are written without coding 686 * code is output_nbits, moved into this routine to reduce overheads 687 * This code potentially breaks if FS>24, so I am limiting 688 * FS to 24 by choice of FSMAX above. 689 */ 690 if (fs > 0) { 691 lbitbuffer <<= fs; 692 lbitbuffer |= v & fsmask; 693 lbits_to_go -= fs; 694 while (lbits_to_go <= 0) { 695 putcbuf((lbitbuffer>>(-lbits_to_go)) & 0xff,buffer); 696 lbits_to_go += 8; 697 } 698 } 699 } 700 /* check if overflowed output buffer */ 701 if (buffer->current > buffer->end) { 702 ffpmsg("rice_encode: end of buffer"); 703 free(diff); 704 return(-1); 232 705 } 233 706 buffer->bitbuffer = lbitbuffer; … … 270 743 int lbitbuffer; 271 744 int lbits_to_go; 745 /* AND mask for the right-most n bits */ 746 static unsigned int mask[33] = 747 {0, 748 0x1, 0x3, 0x7, 0xf, 0x1f, 0x3f, 0x7f, 0xff, 749 0x1ff, 0x3ff, 0x7ff, 0xfff, 0x1fff, 0x3fff, 0x7fff, 0xffff, 750 0x1ffff, 0x3ffff, 0x7ffff, 0xfffff, 0x1fffff, 0x3fffff, 0x7fffff, 0xffffff, 751 0x1ffffff, 0x3ffffff, 0x7ffffff, 0xfffffff, 0x1fffffff, 0x3fffffff, 0x7fffffff, 0xffffffff}; 272 752 273 753 /* … … 282 762 */ 283 763 lbitbuffer <<= lbits_to_go; 284 lbitbuffer |= (bits>>(n-lbits_to_go)) & ((1<<lbits_to_go)-1); 285 if (putcbuf(lbitbuffer & 0xff,buffer) == EOF) return(EOF); 764 /* lbitbuffer |= (bits>>(n-lbits_to_go)) & ((1<<lbits_to_go)-1); */ 765 lbitbuffer |= (bits>>(n-lbits_to_go)) & *(mask+lbits_to_go); 766 putcbuf(lbitbuffer & 0xff,buffer); 286 767 n -= lbits_to_go; 287 768 lbits_to_go = 8; 288 769 } 289 770 lbitbuffer <<= n; 290 lbitbuffer |= ( bits & ((1<<n)-1) ); 771 /* lbitbuffer |= ( bits & ((1<<n)-1) ); */ 772 lbitbuffer |= ( bits & *(mask+n) ); 291 773 lbits_to_go -= n; 292 774 while (lbits_to_go <= 0) { … … 294 776 * bitbuffer full, put out top 8 bits 295 777 */ 296 if (putcbuf((lbitbuffer>>(-lbits_to_go)) & 0xff,buffer) == EOF) 297 return(EOF); 778 putcbuf((lbitbuffer>>(-lbits_to_go)) & 0xff,buffer); 298 779 lbits_to_go += 8; 299 780 } … … 302 783 return(0); 303 784 } 304 305 785 /*---------------------------------------------------------------------------*/ 306 786 /* Flush out the last bits */ … … 309 789 { 310 790 if(buffer->bits_to_go < 8) { 311 if (putcbuf(buffer->bitbuffer<<buffer->bits_to_go,buffer) == EOF) 791 putcbuf(buffer->bitbuffer<<buffer->bits_to_go,buffer); 792 793 /* if (putcbuf(buffer->bitbuffer<<buffer->bits_to_go,buffer) == EOF) 312 794 return(EOF); 795 */ 313 796 } 314 797 return(0); … … 334 817 */ 335 818 819 /*---------------------------------------------------------------------------*/ 336 820 /* this routine used to be called 'rdecomp' (WDP) */ 337 821 … … 356 840 */ 357 841 bsize = 4; 358 /* nblock = 32; */ 842 843 /* nblock = 32; now an input parameter */ 359 844 /* 360 845 * From bsize derive: … … 363 848 * BBITS = bits/pixel for direct coding 364 849 */ 850 851 /* 365 852 switch (bsize) { 366 853 case 1: … … 377 864 break; 378 865 default: 379 f printf (stderr,"rdecomp: bsize must be 1, 2, or 4 bytes");866 ffpmsg("rdecomp: bsize must be 1, 2, or 4 bytes"); 380 867 return 1; 381 868 } 869 */ 870 871 /* move out of switch block, to tweak performance */ 872 fsbits = 5; 873 fsmax = 25; 874 382 875 bbits = 1<<fsbits; 383 876 384 // XXX EAM : sensible code would put this in an init function,385 // with a provided cleanup function. this is just generating386 // the number of 1 bits for any 8bit value387 877 if (nonzero_count == (int *) NULL) { 388 878 /* … … 394 884 nonzero_count = (int *) malloc(256*sizeof(int)); 395 885 if (nonzero_count == (int *) NULL) { 396 f printf (stderr,"rdecomp: insufficient memory");886 ffpmsg("rdecomp: insufficient memory"); 397 887 return 1; 398 888 } … … 411 901 /* first 4 bytes of input buffer contain the value of the first */ 412 902 /* 4 byte integer value, without any encoding */ 413 414 // XXX EAM : wow! this is just doing { lastpix = *(int *)c; } without a cast (encoded value is big-endian)415 903 416 904 lastpix = 0; … … 437 925 } 438 926 fs = (b >> nbits) - 1; 927 439 928 b &= (1<<nbits)-1; 440 929 /* loop over the next block */ … … 494 983 diff = (nzero<<fs) | (b>>nbits); 495 984 b &= (1<<nbits)-1; 985 496 986 /* undo mapping and differencing */ 497 987 if ((diff & 1) == 0) { … … 505 995 } 506 996 if (c > cend) { 507 f printf (stderr,"decompression error: hit end of compressed byte stream");997 ffpmsg("decompression error: hit end of compressed byte stream"); 508 998 return 1; 509 999 } 510 1000 } 511 1001 if (c < cend) { 512 f printf (stderr,"decompression warning: unused bytes at end of compressed buffer");1002 ffpmsg("decompression warning: unused bytes at end of compressed buffer"); 513 1003 } 514 1004 return 0; 515 1005 } 1006 /*---------------------------------------------------------------------------*/ 1007 /* this routine used to be called 'rdecomp' (WDP) */ 1008 1009 int fits_rdecomp_short (unsigned char *c, /* input buffer */ 1010 int clen, /* length of input */ 1011 unsigned short array[], /* output array */ 1012 int nx, /* number of output pixels */ 1013 int nblock) /* coding block size */ 1014 { 1015 int i, imax; 1016 int bsize, k; 1017 int nbits, nzero, fs; 1018 unsigned char *cend, bytevalue; 1019 unsigned int b, diff, lastpix; 1020 int fsmax, fsbits, bbits; 1021 static int *nonzero_count = (int *)NULL; 1022 1023 /* 1024 * Original size of each pixel (bsize, bytes) and coding block 1025 * size (nblock, pixels) 1026 * Could make bsize a parameter to allow more efficient 1027 * compression of short & byte images. 1028 */ 1029 1030 bsize = 2; 1031 1032 /* nblock = 32; now an input parameter */ 1033 /* 1034 * From bsize derive: 1035 * FSBITS = # bits required to store FS 1036 * FSMAX = maximum value for FS 1037 * BBITS = bits/pixel for direct coding 1038 */ 1039 1040 /* 1041 switch (bsize) { 1042 case 1: 1043 fsbits = 3; 1044 fsmax = 6; 1045 break; 1046 case 2: 1047 fsbits = 4; 1048 fsmax = 14; 1049 break; 1050 case 4: 1051 fsbits = 5; 1052 fsmax = 25; 1053 break; 1054 default: 1055 ffpmsg("rdecomp: bsize must be 1, 2, or 4 bytes"); 1056 return 1; 1057 } 1058 */ 1059 1060 /* move out of switch block, to tweak performance */ 1061 fsbits = 4; 1062 fsmax = 14; 1063 1064 bbits = 1<<fsbits; 1065 1066 if (nonzero_count == (int *) NULL) { 1067 /* 1068 * nonzero_count is lookup table giving number of bits 1069 * in 8-bit values not including leading zeros 1070 */ 1071 1072 /* NOTE!!! This memory never gets freed */ 1073 nonzero_count = (int *) malloc(256*sizeof(int)); 1074 if (nonzero_count == (int *) NULL) { 1075 ffpmsg("rdecomp: insufficient memory"); 1076 return 1; 1077 } 1078 nzero = 8; 1079 k = 128; 1080 for (i=255; i>=0; ) { 1081 for ( ; i>=k; i--) nonzero_count[i] = nzero; 1082 k = k/2; 1083 nzero--; 1084 } 1085 } 1086 /* 1087 * Decode in blocks of nblock pixels 1088 */ 1089 1090 /* first 2 bytes of input buffer contain the value of the first */ 1091 /* 2 byte integer value, without any encoding */ 1092 1093 lastpix = 0; 1094 bytevalue = c[0]; 1095 lastpix = lastpix | (bytevalue<<8); 1096 bytevalue = c[1]; 1097 lastpix = lastpix | bytevalue; 1098 1099 c += 2; 1100 cend = c + clen - 2; 1101 1102 b = *c++; /* bit buffer */ 1103 nbits = 8; /* number of bits remaining in b */ 1104 for (i = 0; i<nx; ) { 1105 /* get the FS value from first fsbits */ 1106 nbits -= fsbits; 1107 while (nbits < 0) { 1108 b = (b<<8) | (*c++); 1109 nbits += 8; 1110 } 1111 fs = (b >> nbits) - 1; 1112 1113 b &= (1<<nbits)-1; 1114 /* loop over the next block */ 1115 imax = i + nblock; 1116 if (imax > nx) imax = nx; 1117 if (fs<0) { 1118 /* low-entropy case, all zero differences */ 1119 for ( ; i<imax; i++) array[i] = lastpix; 1120 } else if (fs==fsmax) { 1121 /* high-entropy case, directly coded pixel values */ 1122 for ( ; i<imax; i++) { 1123 k = bbits - nbits; 1124 diff = b<<k; 1125 for (k -= 8; k >= 0; k -= 8) { 1126 b = *c++; 1127 diff |= b<<k; 1128 } 1129 if (nbits>0) { 1130 b = *c++; 1131 diff |= b>>(-k); 1132 b &= (1<<nbits)-1; 1133 } else { 1134 b = 0; 1135 } 1136 1137 /* 1138 * undo mapping and differencing 1139 * Note that some of these operations will overflow the 1140 * unsigned int arithmetic -- that's OK, it all works 1141 * out to give the right answers in the output file. 1142 */ 1143 if ((diff & 1) == 0) { 1144 diff = diff>>1; 1145 } else { 1146 diff = ~(diff>>1); 1147 } 1148 array[i] = diff+lastpix; 1149 lastpix = array[i]; 1150 } 1151 } else { 1152 /* normal case, Rice coding */ 1153 for ( ; i<imax; i++) { 1154 /* count number of leading zeros */ 1155 while (b == 0) { 1156 nbits += 8; 1157 b = *c++; 1158 } 1159 nzero = nbits - nonzero_count[b]; 1160 nbits -= nzero+1; 1161 /* flip the leading one-bit */ 1162 b ^= 1<<nbits; 1163 /* get the FS trailing bits */ 1164 nbits -= fs; 1165 while (nbits < 0) { 1166 b = (b<<8) | (*c++); 1167 nbits += 8; 1168 } 1169 diff = (nzero<<fs) | (b>>nbits); 1170 b &= (1<<nbits)-1; 1171 1172 /* undo mapping and differencing */ 1173 if ((diff & 1) == 0) { 1174 diff = diff>>1; 1175 } else { 1176 diff = ~(diff>>1); 1177 } 1178 array[i] = diff+lastpix; 1179 lastpix = array[i]; 1180 } 1181 } 1182 if (c > cend) { 1183 ffpmsg("decompression error: hit end of compressed byte stream"); 1184 return 1; 1185 } 1186 } 1187 if (c < cend) { 1188 ffpmsg("decompression warning: unused bytes at end of compressed buffer"); 1189 } 1190 return 0; 1191 } 1192 /*---------------------------------------------------------------------------*/ 1193 /* this routine used to be called 'rdecomp' (WDP) */ 1194 1195 int fits_rdecomp_byte (unsigned char *c, /* input buffer */ 1196 int clen, /* length of input */ 1197 unsigned char array[], /* output array */ 1198 int nx, /* number of output pixels */ 1199 int nblock) /* coding block size */ 1200 { 1201 int i, imax; 1202 int bsize, k; 1203 int nbits, nzero, fs; 1204 unsigned char *cend; 1205 unsigned int b, diff, lastpix; 1206 int fsmax, fsbits, bbits; 1207 static int *nonzero_count = (int *)NULL; 1208 1209 /* 1210 * Original size of each pixel (bsize, bytes) and coding block 1211 * size (nblock, pixels) 1212 * Could make bsize a parameter to allow more efficient 1213 * compression of short & byte images. 1214 */ 1215 1216 bsize = 1; 1217 1218 /* nblock = 32; now an input parameter */ 1219 /* 1220 * From bsize derive: 1221 * FSBITS = # bits required to store FS 1222 * FSMAX = maximum value for FS 1223 * BBITS = bits/pixel for direct coding 1224 */ 1225 1226 /* 1227 switch (bsize) { 1228 case 1: 1229 fsbits = 3; 1230 fsmax = 6; 1231 break; 1232 case 2: 1233 fsbits = 4; 1234 fsmax = 14; 1235 break; 1236 case 4: 1237 fsbits = 5; 1238 fsmax = 25; 1239 break; 1240 default: 1241 ffpmsg("rdecomp: bsize must be 1, 2, or 4 bytes"); 1242 return 1; 1243 } 1244 */ 1245 1246 /* move out of switch block, to tweak performance */ 1247 fsbits = 3; 1248 fsmax = 6; 1249 1250 bbits = 1<<fsbits; 1251 1252 if (nonzero_count == (int *) NULL) { 1253 /* 1254 * nonzero_count is lookup table giving number of bits 1255 * in 8-bit values not including leading zeros 1256 */ 1257 1258 /* NOTE!!! This memory never gets freed */ 1259 nonzero_count = (int *) malloc(256*sizeof(int)); 1260 if (nonzero_count == (int *) NULL) { 1261 ffpmsg("rdecomp: insufficient memory"); 1262 return 1; 1263 } 1264 nzero = 8; 1265 k = 128; 1266 for (i=255; i>=0; ) { 1267 for ( ; i>=k; i--) nonzero_count[i] = nzero; 1268 k = k/2; 1269 nzero--; 1270 } 1271 } 1272 /* 1273 * Decode in blocks of nblock pixels 1274 */ 1275 1276 /* first byte of input buffer contain the value of the first */ 1277 /* byte integer value, without any encoding */ 1278 1279 lastpix = c[0]; 1280 c += 1; 1281 cend = c + clen - 1; 1282 1283 b = *c++; /* bit buffer */ 1284 nbits = 8; /* number of bits remaining in b */ 1285 for (i = 0; i<nx; ) { 1286 /* get the FS value from first fsbits */ 1287 nbits -= fsbits; 1288 while (nbits < 0) { 1289 b = (b<<8) | (*c++); 1290 nbits += 8; 1291 } 1292 fs = (b >> nbits) - 1; 1293 1294 b &= (1<<nbits)-1; 1295 /* loop over the next block */ 1296 imax = i + nblock; 1297 if (imax > nx) imax = nx; 1298 if (fs<0) { 1299 /* low-entropy case, all zero differences */ 1300 for ( ; i<imax; i++) array[i] = lastpix; 1301 } else if (fs==fsmax) { 1302 /* high-entropy case, directly coded pixel values */ 1303 for ( ; i<imax; i++) { 1304 k = bbits - nbits; 1305 diff = b<<k; 1306 for (k -= 8; k >= 0; k -= 8) { 1307 b = *c++; 1308 diff |= b<<k; 1309 } 1310 if (nbits>0) { 1311 b = *c++; 1312 diff |= b>>(-k); 1313 b &= (1<<nbits)-1; 1314 } else { 1315 b = 0; 1316 } 1317 1318 /* 1319 * undo mapping and differencing 1320 * Note that some of these operations will overflow the 1321 * unsigned int arithmetic -- that's OK, it all works 1322 * out to give the right answers in the output file. 1323 */ 1324 if ((diff & 1) == 0) { 1325 diff = diff>>1; 1326 } else { 1327 diff = ~(diff>>1); 1328 } 1329 array[i] = diff+lastpix; 1330 lastpix = array[i]; 1331 } 1332 } else { 1333 /* normal case, Rice coding */ 1334 for ( ; i<imax; i++) { 1335 /* count number of leading zeros */ 1336 while (b == 0) { 1337 nbits += 8; 1338 b = *c++; 1339 } 1340 nzero = nbits - nonzero_count[b]; 1341 nbits -= nzero+1; 1342 /* flip the leading one-bit */ 1343 b ^= 1<<nbits; 1344 /* get the FS trailing bits */ 1345 nbits -= fs; 1346 while (nbits < 0) { 1347 b = (b<<8) | (*c++); 1348 nbits += 8; 1349 } 1350 diff = (nzero<<fs) | (b>>nbits); 1351 b &= (1<<nbits)-1; 1352 1353 /* undo mapping and differencing */ 1354 if ((diff & 1) == 0) { 1355 diff = diff>>1; 1356 } else { 1357 diff = ~(diff>>1); 1358 } 1359 array[i] = diff+lastpix; 1360 lastpix = array[i]; 1361 } 1362 } 1363 if (c > cend) { 1364 ffpmsg("decompression error: hit end of compressed byte stream"); 1365 return 1; 1366 } 1367 } 1368 if (c < cend) { 1369 ffpmsg("decompression warning: unused bytes at end of compressed buffer"); 1370 } 1371 return 0; 1372 }
Note:
See TracChangeset
for help on using the changeset viewer.
