IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 18269


Ignore:
Timestamp:
Jun 22, 2008, 8:52:54 AM (18 years ago)
Author:
eugene
Message:

replaced with CFITSIO updates from 3.080; defined our own ricecomp.h to make upgrades easier

Location:
trunk/Ohana/src/libfits/extern
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/libfits/extern/ricecomp.c

    r15487 r18269  
     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
    17/*
    28  The following code was written by Richard White at STScI and made
    … …  
    2329#include <stdlib.h>
    2430#include <string.h>
    25 #include "ricecomp.h"  /* originally included in rcomp.c file (WDP) */
     31
     32typedef unsigned char Buffer_t;
     33
     34typedef 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
    2644// #include "fitsio2.h"
    27 
     45# include "ricecomp.h"
    2846
    2947static void start_outputing_bits(Buffer *buffer);
    … …  
    3250
    3351/* this routine used to be called 'rcomp'  (WDP)  */
     52/*---------------------------------------------------------------------------*/
    3453
    3554int fits_rcomp(int a[],         /* input array                  */
    … …  
    5574     */
    5675    bsize = 4;
    57 /*    nblock = 32; */
     76
     77/*    nblock = 32; now an input parameter*/
    5878    /*
    5979     * From bsize derive:
    … …  
    6282     * BBITS = bits/pixel for direct coding
    6383     */
     84
     85/*
    6486    switch (bsize) {
    6587    case 1:
    … …  
    7698        break;
    7799    default:
    78         fprintf (stderr, "rdecomp: bsize must be 1, 2, or 4 bytes");
     100        ffpmsg("rdecomp: bsize must be 1, 2, or 4 bytes");
    79101        return(-1);
    80102    }
     103*/
     104
     105    /* move out of switch block, to tweak performance */
     106    fsbits = 5;
     107    fsmax = 25;
    81108    bbits = 1<<fsbits;
    82109
    … …  
    93120    diff = (unsigned int *) malloc(nblock*sizeof(unsigned int));
    94121    if (diff == (unsigned int *) NULL) {
    95         fprintf (stderr, "fits_rcomp: insufficient memory");
     122        ffpmsg("fits_rcomp: insufficient memory");
    96123        return(-1);
    97124    }
    … …  
    103130    /* write out first int value to the first 4 bytes of the buffer */
    104131    if (output_nbits(buffer, a[0], 32) == EOF) {
    105         fprintf (stderr, "rice_encode: end of buffer");
     132        ffpmsg("rice_encode: end of buffer");
    106133        free(diff);
    107134        return(-1);
    … …  
    133160            lastpix = nextpix;
    134161        }
     162
    135163        /*
    136164         * compute number of bits to split from sum
    … …  
    140168        psum = ((unsigned int) dpsum ) >> 1;
    141169        for (fs = 0; psum>0; fs++) psum >>= 1;
     170
    142171        /*
    143172         * write the codes
    … …  
    149178             */
    150179            if (output_nbits(buffer, fsmax+1, fsbits) == EOF) {
    151                 fprintf (stderr, "rice_encode: end of buffer");
     180                ffpmsg("rice_encode: end of buffer");
    152181                free(diff);
    153182                return(-1);
    … …  
    155184            for (j=0; j<thisblock; j++) {
    156185                if (output_nbits(buffer, diff[j], bbits) == EOF) {
    157                     fprintf (stderr, "rice_encode: end of buffer");
     186                    ffpmsg("rice_encode: end of buffer");
    158187                    free(diff);
    159188                    return(-1);
    … …  
    167196             */
    168197            if (output_nbits(buffer, 0, fsbits) == EOF) {
    169                 fprintf (stderr, "rice_encode: end of buffer");
     198                ffpmsg("rice_encode: end of buffer");
    170199                free(diff);
    171200                return(-1);
    … …  
    174203            /* normal case: not either very high or very low entropy */
    175204            if (output_nbits(buffer, fs+1, fsbits) == EOF) {
    176                 fprintf (stderr, "rice_encode: end of buffer");
     205                ffpmsg("rice_encode: end of buffer");
    177206                free(diff);
    178207                return(-1);
    … …  
    196225                } else {
    197226                    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
    203229                    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);
    209231                    }
    210232                    lbitbuffer = 1;
    … …  
    222244                    lbits_to_go -= fs;
    223245                    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);
    229247                        lbits_to_go += 8;
    230248                    }
    231249                }
     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
     271int 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{
     278Buffer bufmem, *buffer = &bufmem;
     279int bsize, i, j, thisblock;
     280
     281/*
     282NOTE: in principle, the following 2 variable could be declared as 'short'
     283but in fact the code runs faster (on 32-bit Linux at least) as 'int'
     284*/
     285int lastpix, nextpix;
     286/* int pdiff; */
     287short pdiff;
     288int v, fs, fsmask, top, fsmax, fsbits, bbits;
     289int lbitbuffer, lbits_to_go;
     290/* unsigned int psum; */
     291unsigned short psum;
     292double pixelsum, dpsum;
     293unsigned 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
     495int 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{
     502Buffer bufmem, *buffer = &bufmem;
     503int bsize, i, j, thisblock;
     504
     505/*
     506NOTE: in principle, the following 2 variable could be declared as 'short'
     507but in fact the code runs faster (on 32-bit Linux at least) as 'int'
     508*/
     509int lastpix, nextpix;
     510/* int pdiff; */
     511signed char pdiff;
     512int v, fs, fsmask, top, fsmax, fsbits, bbits;
     513int lbitbuffer, lbits_to_go;
     514/* unsigned int psum; */
     515unsigned char psum;
     516double pixelsum, dpsum;
     517unsigned 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);
    232705            }
    233706            buffer->bitbuffer = lbitbuffer;
    … …  
    270743int lbitbuffer;
    271744int 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};
    272752
    273753    /*
    … …  
    282762         */
    283763        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);
    286767        n -= lbits_to_go;
    287768        lbits_to_go = 8;
    288769    }
    289770    lbitbuffer <<= n;
    290     lbitbuffer |= ( bits & ((1<<n)-1) );
     771/*    lbitbuffer |= ( bits & ((1<<n)-1) ); */
     772    lbitbuffer |= ( bits & *(mask+n) );
    291773    lbits_to_go -= n;
    292774    while (lbits_to_go <= 0) {
    … …  
    294776         * bitbuffer full, put out top 8 bits
    295777         */
    296         if (putcbuf((lbitbuffer>>(-lbits_to_go)) & 0xff,buffer) == EOF)
    297             return(EOF);
     778        putcbuf((lbitbuffer>>(-lbits_to_go)) & 0xff,buffer);
    298779        lbits_to_go += 8;
    299780    }
    … …  
    302783    return(0);
    303784}
    304 
    305785/*---------------------------------------------------------------------------*/
    306786/* Flush out the last bits */
    … …  
    309789{
    310790    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)
    312794            return(EOF);
     795*/
    313796    }
    314797    return(0);
    … …  
    334817*/
    335818
     819/*---------------------------------------------------------------------------*/
    336820/* this routine used to be called 'rdecomp'  (WDP)  */
    337821
    … …  
    356840     */
    357841    bsize = 4;
    358 /*    nblock = 32; */
     842
     843/*    nblock = 32; now an input parameter */
    359844    /*
    360845     * From bsize derive:
    … …  
    363848     * BBITS = bits/pixel for direct coding
    364849     */
     850
     851/*
    365852    switch (bsize) {
    366853    case 1:
    … …  
    377864        break;
    378865    default:
    379         fprintf (stderr, "rdecomp: bsize must be 1, 2, or 4 bytes");
     866        ffpmsg("rdecomp: bsize must be 1, 2, or 4 bytes");
    380867        return 1;
    381868    }
     869*/
     870
     871    /* move out of switch block, to tweak performance */
     872    fsbits = 5;
     873    fsmax = 25;
     874
    382875    bbits = 1<<fsbits;
    383876
    384     // XXX EAM : sensible code would put this in an init function,
    385     // with a provided cleanup function.  this is just generating
    386     // the number of 1 bits for any 8bit value
    387877    if (nonzero_count == (int *) NULL) {
    388878        /*
    … …  
    394884        nonzero_count = (int *) malloc(256*sizeof(int));
    395885        if (nonzero_count == (int *) NULL) {
    396             fprintf (stderr, "rdecomp: insufficient memory");
     886            ffpmsg("rdecomp: insufficient memory");
    397887            return 1;
    398888        }
    … …  
    411901    /* first 4 bytes of input buffer contain the value of the first */
    412902    /* 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)
    415903   
    416904    lastpix = 0;
    … …  
    437925        }
    438926        fs = (b >> nbits) - 1;
     927
    439928        b &= (1<<nbits)-1;
    440929        /* loop over the next block */
    … …  
    494983                diff = (nzero<<fs) | (b>>nbits);
    495984                b &= (1<<nbits)-1;
     985
    496986                /* undo mapping and differencing */
    497987                if ((diff & 1) == 0) {
    … …  
    505995        }
    506996        if (c > cend) {
    507             fprintf (stderr, "decompression error: hit end of compressed byte stream");
     997            ffpmsg("decompression error: hit end of compressed byte stream");
    508998            return 1;
    509999        }
    5101000    }
    5111001    if (c < cend) {
    512         fprintf (stderr, "decompression warning: unused bytes at end of compressed buffer");
     1002        ffpmsg("decompression warning: unused bytes at end of compressed buffer");
    5131003    }
    5141004    return 0;
    5151005}
     1006/*---------------------------------------------------------------------------*/
     1007/* this routine used to be called 'rdecomp'  (WDP)  */
     1008
     1009int 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{
     1015int i, imax;
     1016int bsize, k;
     1017int nbits, nzero, fs;
     1018unsigned char *cend, bytevalue;
     1019unsigned int b, diff, lastpix;
     1020int fsmax, fsbits, bbits;
     1021static 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
     1195int 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{
     1201int i, imax;
     1202int bsize, k;
     1203int nbits, nzero, fs;
     1204unsigned char *cend;
     1205unsigned int b, diff, lastpix;
     1206int fsmax, fsbits, bbits;
     1207static 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}

  • trunk/Ohana/src/libfits/extern/ricecomp.h

    r15487 r18269  
    1 /* @(#) buffer.h 1.1 98/07/21 12:34:27 */
    2 /* buffer.h: structure for compression to buffer rather than to a file, including
    3  *   bit I/O buffer
    4  *
    5  * R. White, 19 June 1998
    6  */
    71
    82
    9 typedef unsigned char Buffer_t;
     3// 2008.06.22 EAM : This code was taken from CFITSIO and included in Ohana for rice
     4// decompression.  In order to include this .c file (or future upgrades), we include our own
     5// version of ricecomp.h (not the one available in the CFITSIO tree) and include it instead of
     6// fitsio2.h.  This file defines our own version of ffpmsg. 
    107
    11 typedef struct {
    12         int bitbuffer;          /* bit buffer                                   */
    13         int bits_to_go;         /* bits to go in buffer                 */
    14         Buffer_t *start;        /* start of buffer                              */
    15         Buffer_t *current;      /* current position in buffer   */
    16         Buffer_t *end;          /* end of buffer                                */
    17 } Buffer;
    18 
    19 #define buffree(mf)             (free(mf->start), free(mf))
    20 #define bufused(mf)             (mf->current - mf->start)
    21 #define bufreset(mf)    (mf->current = mf->start)
    22 
    23 /*
    24  * getcbuf, putcbuf macros for character IO to buffer
    25  * putcbuf returns EOF on end of buffer, else returns 0
    26  */
    27 #define getcbuf(mf) ((mf->current >= mf->end) ? EOF : *(mf->current)++)
    28 #define putcbuf(c,mf) \
    29         ((mf->current >= mf->end) ? \
    30         EOF :\
    31         ((*(mf->current)++ = c), 0))
    32 
    33 /*
    34  * bufalloc sets up buffer of length n
    35  */
    36 
    37 /*  not needed by CFITSIO
    38 
    39 static Buffer *bufalloc(int n)
    40 {
    41 Buffer *mf;
    42 
    43         mf = (Buffer *) malloc(sizeof(Buffer));
    44         if (mf == (Buffer *)NULL) return((Buffer *)NULL);
    45 
    46         mf->start = (Buffer_t *) malloc(n*sizeof(Buffer_t));
    47         if (mf->start == (Buffer_t *)NULL) {
    48                 free(mf);
    49                 return((Buffer *)NULL);
    50         }
    51         mf->bits_to_go = 8;
    52         mf->end = mf->start + n;
    53         mf->current = mf->start;
    54         return(mf);
     8int ffpmsg (char *line) {
     9   
     10    fprintf (stderr, "%s\n", line);
     11    return 1;
    5512}
    56 */
    57 
    58 /*
    59  * bufrealloc extends buffer (or truncates it) by
    60  * reallocating memory
    61  */
    62 
    63 /* not needed by CFITSIO
    64 static int bufrealloc(Buffer *mf, int n)
    65 {
    66 int len;
    67 
    68         len = mf->current - mf->start;
    69 
    70         * silently throw away data if buffer is already longer than n *
    71         if (len>n) len = n;
    72         if (len<0) len = 0;
    73 
    74         mf->start = (Buffer_t *) realloc(mf->start, n*sizeof(Buffer_t));
    75         if (mf->start == (Buffer_t *)NULL) return(0);
    76 
    77         mf->end = mf->start + n;
    78         mf->current = mf->start + len;
    79         return(n);
    80 }
    81 */
    82 
    83 /*
    84  * bufdump dumps contents of buffer to outfile and resets
    85  * it to be empty.  Returns number of bytes written.
    86  *
    87  * Note we don't write out the bit buffer -- you must call
    88  * done_outputing_bits() first to ensure that the bit buffer
    89  * is written out.  I do it this way to allow incremental
    90  * buffer dumps while bit IO is still going on.
    91  */
    92 
    93 /*  not needed by CFITSIO
    94 
    95 static int bufdump(FILE *outfile, Buffer *buffer)
    96 {
    97 int ndump;
    98 
    99         ndump = bufused(buffer);
    100         if (fwrite(buffer->start, 1, ndump, outfile) != ndump) {
    101                 fprintf(stderr, "bufdump: error in write\n");
    102                 exit(1);
    103     }
    104         bufreset(buffer);
    105         return(ndump);
    106 }
    107 */
Note: See TracChangeset for help on using the changeset viewer.