Index: /trunk/Ohana/src/libfits/extern/ricecomp.c
===================================================================
--- /trunk/Ohana/src/libfits/extern/ricecomp.c	(revision 18268)
+++ /trunk/Ohana/src/libfits/extern/ricecomp.c	(revision 18269)
@@ -1,2 +1,8 @@
+
+// 2008.06.22 EAM : This code was taken from CFITSIO and included in Ohana for rice
+// decompression.  In order to include this .c file (or future upgrades), we include our own
+// version of ricecomp.h (not the one available in the CFITSIO tree) and include it instead of
+// fitsio2.h.  This file defines our own version of ffpmsg.  
+
 /*
   The following code was written by Richard White at STScI and made
@@ -23,7 +29,19 @@
 #include <stdlib.h>
 #include <string.h>
-#include "ricecomp.h"  /* originally included in rcomp.c file (WDP) */
+
+typedef unsigned char Buffer_t;
+
+typedef struct {
+	int bitbuffer;		/* bit buffer					*/
+	int bits_to_go;		/* bits to go in buffer			*/
+	Buffer_t *start;	/* start of buffer				*/
+	Buffer_t *current;	/* current position in buffer	*/
+	Buffer_t *end;		/* end of buffer				*/
+} Buffer;
+
+#define putcbuf(c,mf) 	((*(mf->current)++ = c), 0)
+
 // #include "fitsio2.h"
-
+# include "ricecomp.h"
 
 static void start_outputing_bits(Buffer *buffer);
@@ -32,4 +50,5 @@
 
 /* this routine used to be called 'rcomp'  (WDP)  */
+/*---------------------------------------------------------------------------*/
 
 int fits_rcomp(int a[],		/* input array			*/
@@ -55,5 +74,6 @@
      */
     bsize = 4;
-/*    nblock = 32; */
+
+/*    nblock = 32; now an input parameter*/
     /*
      * From bsize derive:
@@ -62,4 +82,6 @@
      * BBITS = bits/pixel for direct coding
      */
+
+/*
     switch (bsize) {
     case 1:
@@ -76,7 +98,12 @@
 	break;
     default:
-        fprintf (stderr, "rdecomp: bsize must be 1, 2, or 4 bytes");
+        ffpmsg("rdecomp: bsize must be 1, 2, or 4 bytes");
 	return(-1);
     }
+*/
+
+    /* move out of switch block, to tweak performance */
+    fsbits = 5;
+    fsmax = 25;
     bbits = 1<<fsbits;
 
@@ -93,5 +120,5 @@
     diff = (unsigned int *) malloc(nblock*sizeof(unsigned int));
     if (diff == (unsigned int *) NULL) {
-        fprintf (stderr, "fits_rcomp: insufficient memory");
+        ffpmsg("fits_rcomp: insufficient memory");
 	return(-1);
     }
@@ -103,5 +130,5 @@
     /* write out first int value to the first 4 bytes of the buffer */
     if (output_nbits(buffer, a[0], 32) == EOF) {
-        fprintf (stderr, "rice_encode: end of buffer");
+        ffpmsg("rice_encode: end of buffer");
         free(diff);
         return(-1);
@@ -133,4 +160,5 @@
 	    lastpix = nextpix;
 	}
+
 	/*
 	 * compute number of bits to split from sum
@@ -140,4 +168,5 @@
 	psum = ((unsigned int) dpsum ) >> 1;
 	for (fs = 0; psum>0; fs++) psum >>= 1;
+
 	/*
 	 * write the codes
@@ -149,5 +178,5 @@
 	     */
 	    if (output_nbits(buffer, fsmax+1, fsbits) == EOF) {
-                fprintf (stderr, "rice_encode: end of buffer");
+                ffpmsg("rice_encode: end of buffer");
                 free(diff);
 		return(-1);
@@ -155,5 +184,5 @@
 	    for (j=0; j<thisblock; j++) {
 		if (output_nbits(buffer, diff[j], bbits) == EOF) {
-                    fprintf (stderr, "rice_encode: end of buffer");
+                    ffpmsg("rice_encode: end of buffer");
                     free(diff);
 		    return(-1);
@@ -167,5 +196,5 @@
 	     */
 	    if (output_nbits(buffer, 0, fsbits) == EOF) {
-                fprintf (stderr, "rice_encode: end of buffer");
+                ffpmsg("rice_encode: end of buffer");
                 free(diff);
 		return(-1);
@@ -174,5 +203,5 @@
 	    /* normal case: not either very high or very low entropy */
 	    if (output_nbits(buffer, fs+1, fsbits) == EOF) {
-                fprintf (stderr, "rice_encode: end of buffer");
+                ffpmsg("rice_encode: end of buffer");
                 free(diff);
 		return(-1);
@@ -196,15 +225,8 @@
 		} else {
 		    lbitbuffer <<= lbits_to_go;
-		    if (putcbuf(lbitbuffer & 0xff,buffer) == EOF) {
-                        fprintf (stderr, "rice_encode: end of buffer");
-                        free(diff);
-			return(-1);
-		    }
+		    putcbuf(lbitbuffer & 0xff,buffer);
+
 		    for (top -= lbits_to_go; top>=8; top -= 8) {
-			if (putcbuf(0, buffer) == EOF) {
-                            fprintf (stderr, "rice_encode: end of buffer");
-                            free(diff);
-			    return(-1);
-			}
+			putcbuf(0, buffer);
 		    }
 		    lbitbuffer = 1;
@@ -222,12 +244,463 @@
 		    lbits_to_go -= fs;
 		    while (lbits_to_go <= 0) {
-			if (putcbuf((lbitbuffer>>(-lbits_to_go)) & 0xff,buffer)==EOF) {
-                            fprintf (stderr, "rice_encode: end of buffer");
-                            free(diff);
-			    return(-1);
-			}
+			putcbuf((lbitbuffer>>(-lbits_to_go)) & 0xff,buffer);
 			lbits_to_go += 8;
 		    }
 		}
+	    }
+
+	    /* check if overflowed output buffer */
+	    if (buffer->current > buffer->end) {
+                 ffpmsg("rice_encode: end of buffer");
+                 free(diff);
+		 return(-1);
+	    }
+	    buffer->bitbuffer = lbitbuffer;
+	    buffer->bits_to_go = lbits_to_go;
+	}
+    }
+    done_outputing_bits(buffer);
+    free(diff);
+    /*
+     * return number of bytes used
+     */
+    return(buffer->current - buffer->start);
+}
+/*---------------------------------------------------------------------------*/
+
+int fits_rcomp_short(
+	  short a[],		/* input array			*/
+	  int nx,		/* number of input pixels	*/
+	  unsigned char *c,	/* output buffer		*/
+	  int clen,		/* max length of output		*/
+	  int nblock)		/* coding block size		*/
+{
+Buffer bufmem, *buffer = &bufmem;
+int bsize, i, j, thisblock;
+
+/* 
+NOTE: in principle, the following 2 variable could be declared as 'short'
+but in fact the code runs faster (on 32-bit Linux at least) as 'int'
+*/
+int lastpix, nextpix;
+/* int pdiff; */
+short pdiff; 
+int v, fs, fsmask, top, fsmax, fsbits, bbits;
+int lbitbuffer, lbits_to_go;
+/* unsigned int psum; */
+unsigned short psum;
+double pixelsum, dpsum;
+unsigned int *diff;
+
+    /*
+     * Original size of each pixel (bsize, bytes) and coding block
+     * size (nblock, pixels)
+     * Could make bsize a parameter to allow more efficient
+     * compression of short & byte images.
+     */
+    bsize = 2;
+
+/*    nblock = 32; now an input parameter */
+    /*
+     * From bsize derive:
+     * FSBITS = # bits required to store FS
+     * FSMAX = maximum value for FS
+     * BBITS = bits/pixel for direct coding
+     */
+
+/*
+    switch (bsize) {
+    case 1:
+	fsbits = 3;
+	fsmax = 6;
+	break;
+    case 2:
+	fsbits = 4;
+	fsmax = 14;
+	break;
+    case 4:
+	fsbits = 5;
+	fsmax = 25;
+	break;
+    default:
+        ffpmsg("rdecomp: bsize must be 1, 2, or 4 bytes");
+	return(-1);
+    }
+*/
+
+    /* move these out of switch block to further tweak performance */
+    fsbits = 4;
+    fsmax = 14;
+    bbits = 1<<fsbits;
+
+    /*
+     * Set up buffer pointers
+     */
+    buffer->start = c;
+    buffer->current = c;
+    buffer->end = c+clen;
+    buffer->bits_to_go = 8;
+    /*
+     * array for differences mapped to non-negative values
+     */
+    diff = (unsigned int *) malloc(nblock*sizeof(unsigned int));
+    if (diff == (unsigned int *) NULL) {
+        ffpmsg("fits_rcomp: insufficient memory");
+	return(-1);
+    }
+    /*
+     * Code in blocks of nblock pixels
+     */
+    start_outputing_bits(buffer);
+
+    /* write out first short value to the first 2 bytes of the buffer */
+    if (output_nbits(buffer, a[0], 16) == EOF) {
+        ffpmsg("rice_encode: end of buffer");
+        free(diff);
+        return(-1);
+    }
+
+    lastpix = a[0];  /* the first difference will always be zero */
+
+    thisblock = nblock;
+    for (i=0; i<nx; i += nblock) {
+	/* last block may be shorter */
+	if (nx-i < nblock) thisblock = nx-i;
+	/*
+	 * Compute differences of adjacent pixels and map them to unsigned values.
+	 * Note that this may overflow the integer variables -- that's
+	 * OK, because we can recover when decompressing.  If we were
+	 * compressing shorts or bytes, would want to do this arithmetic
+	 * with short/byte working variables (though diff will still be
+	 * passed as an int.)
+	 *
+	 * compute sum of mapped pixel values at same time
+	 * use double precision for sum to allow 32-bit integer inputs
+	 */
+	pixelsum = 0.0;
+	for (j=0; j<thisblock; j++) {
+	    nextpix = a[i+j];
+	    pdiff = nextpix - lastpix;
+	    diff[j] = (unsigned int) ((pdiff<0) ? ~(pdiff<<1) : (pdiff<<1));
+	    pixelsum += diff[j];
+	    lastpix = nextpix;
+	}
+	/*
+	 * compute number of bits to split from sum
+	 */
+	dpsum = (pixelsum - (thisblock/2) - 1)/thisblock;
+	if (dpsum < 0) dpsum = 0.0;
+/*	psum = ((unsigned int) dpsum ) >> 1; */
+	psum = ((unsigned short) dpsum ) >> 1;
+	for (fs = 0; psum>0; fs++) psum >>= 1;
+
+	/*
+	 * write the codes
+	 * fsbits ID bits used to indicate split level
+	 */
+	if (fs >= fsmax) {
+	    /* Special high entropy case when FS >= fsmax
+	     * Just write pixel difference values directly, no Rice coding at all.
+	     */
+	    if (output_nbits(buffer, fsmax+1, fsbits) == EOF) {
+                ffpmsg("rice_encode: end of buffer");
+                free(diff);
+		return(-1);
+	    }
+	    for (j=0; j<thisblock; j++) {
+		if (output_nbits(buffer, diff[j], bbits) == EOF) {
+                    ffpmsg("rice_encode: end of buffer");
+                    free(diff);
+		    return(-1);
+		}
+	    }
+	} else if (fs == 0 && pixelsum == 0) {
+	    /*
+	     * special low entropy case when FS = 0 and pixelsum=0 (all
+	     * pixels in block are zero.)
+	     * Output a 0 and return
+	     */
+	    if (output_nbits(buffer, 0, fsbits) == EOF) {
+                ffpmsg("rice_encode: end of buffer");
+                free(diff);
+		return(-1);
+	    }
+	} else {
+	    /* normal case: not either very high or very low entropy */
+	    if (output_nbits(buffer, fs+1, fsbits) == EOF) {
+                ffpmsg("rice_encode: end of buffer");
+                free(diff);
+		return(-1);
+	    }
+	    fsmask = (1<<fs) - 1;
+	    /*
+	     * local copies of bit buffer to improve optimization
+	     */
+	    lbitbuffer = buffer->bitbuffer;
+	    lbits_to_go = buffer->bits_to_go;
+	    for (j=0; j<thisblock; j++) {
+		v = diff[j];
+		top = v >> fs;
+		/*
+		 * top is coded by top zeros + 1
+		 */
+		if (lbits_to_go >= top+1) {
+		    lbitbuffer <<= top+1;
+		    lbitbuffer |= 1;
+		    lbits_to_go -= top+1;
+		} else {
+		    lbitbuffer <<= lbits_to_go;
+		    putcbuf(lbitbuffer & 0xff,buffer);
+		    for (top -= lbits_to_go; top>=8; top -= 8) {
+			putcbuf(0, buffer);
+		    }
+		    lbitbuffer = 1;
+		    lbits_to_go = 7-top;
+		}
+		/*
+		 * bottom FS bits are written without coding
+		 * code is output_nbits, moved into this routine to reduce overheads
+		 * This code potentially breaks if FS>24, so I am limiting
+		 * FS to 24 by choice of FSMAX above.
+		 */
+		if (fs > 0) {
+		    lbitbuffer <<= fs;
+		    lbitbuffer |= v & fsmask;
+		    lbits_to_go -= fs;
+		    while (lbits_to_go <= 0) {
+			putcbuf((lbitbuffer>>(-lbits_to_go)) & 0xff,buffer);
+			lbits_to_go += 8;
+		    }
+		}
+	    }
+	    /* check if overflowed output buffer */
+	    if (buffer->current > buffer->end) {
+                 ffpmsg("rice_encode: end of buffer");
+                 free(diff);
+		 return(-1);
+	    }
+	    buffer->bitbuffer = lbitbuffer;
+	    buffer->bits_to_go = lbits_to_go;
+	}
+    }
+    done_outputing_bits(buffer);
+    free(diff);
+    /*
+     * return number of bytes used
+     */
+    return(buffer->current - buffer->start);
+}
+/*---------------------------------------------------------------------------*/
+
+int fits_rcomp_byte(
+	  signed char a[],		/* input array			*/
+	  int nx,		/* number of input pixels	*/
+	  unsigned char *c,	/* output buffer		*/
+	  int clen,		/* max length of output		*/
+	  int nblock)		/* coding block size		*/
+{
+Buffer bufmem, *buffer = &bufmem;
+int bsize, i, j, thisblock;
+
+/* 
+NOTE: in principle, the following 2 variable could be declared as 'short'
+but in fact the code runs faster (on 32-bit Linux at least) as 'int'
+*/
+int lastpix, nextpix;
+/* int pdiff; */
+signed char pdiff; 
+int v, fs, fsmask, top, fsmax, fsbits, bbits;
+int lbitbuffer, lbits_to_go;
+/* unsigned int psum; */
+unsigned char psum;
+double pixelsum, dpsum;
+unsigned int *diff;
+
+    /*
+     * Original size of each pixel (bsize, bytes) and coding block
+     * size (nblock, pixels)
+     * Could make bsize a parameter to allow more efficient
+     * compression of short & byte images.
+     */
+    bsize = 1;
+
+/*    nblock = 32; now an input parameter */
+    /*
+     * From bsize derive:
+     * FSBITS = # bits required to store FS
+     * FSMAX = maximum value for FS
+     * BBITS = bits/pixel for direct coding
+     */
+
+/*
+    switch (bsize) {
+    case 1:
+	fsbits = 3;
+	fsmax = 6;
+	break;
+    case 2:
+	fsbits = 4;
+	fsmax = 14;
+	break;
+    case 4:
+	fsbits = 5;
+	fsmax = 25;
+	break;
+    default:
+        ffpmsg("rdecomp: bsize must be 1, 2, or 4 bytes");
+	return(-1);
+    }
+*/
+
+    /* move these out of switch block to further tweak performance */
+    fsbits = 3;
+    fsmax = 6;
+    bbits = 1<<fsbits;
+
+    /*
+     * Set up buffer pointers
+     */
+    buffer->start = c;
+    buffer->current = c;
+    buffer->end = c+clen;
+    buffer->bits_to_go = 8;
+    /*
+     * array for differences mapped to non-negative values
+     */
+    diff = (unsigned int *) malloc(nblock*sizeof(unsigned int));
+    if (diff == (unsigned int *) NULL) {
+        ffpmsg("fits_rcomp: insufficient memory");
+	return(-1);
+    }
+    /*
+     * Code in blocks of nblock pixels
+     */
+    start_outputing_bits(buffer);
+
+    /* write out first byte value to the first  byte of the buffer */
+    if (output_nbits(buffer, a[0], 8) == EOF) {
+        ffpmsg("rice_encode: end of buffer");
+        free(diff);
+        return(-1);
+    }
+
+    lastpix = a[0];  /* the first difference will always be zero */
+
+    thisblock = nblock;
+    for (i=0; i<nx; i += nblock) {
+	/* last block may be shorter */
+	if (nx-i < nblock) thisblock = nx-i;
+	/*
+	 * Compute differences of adjacent pixels and map them to unsigned values.
+	 * Note that this may overflow the integer variables -- that's
+	 * OK, because we can recover when decompressing.  If we were
+	 * compressing shorts or bytes, would want to do this arithmetic
+	 * with short/byte working variables (though diff will still be
+	 * passed as an int.)
+	 *
+	 * compute sum of mapped pixel values at same time
+	 * use double precision for sum to allow 32-bit integer inputs
+	 */
+	pixelsum = 0.0;
+	for (j=0; j<thisblock; j++) {
+	    nextpix = a[i+j];
+	    pdiff = nextpix - lastpix;
+	    diff[j] = (unsigned int) ((pdiff<0) ? ~(pdiff<<1) : (pdiff<<1));
+	    pixelsum += diff[j];
+	    lastpix = nextpix;
+	}
+	/*
+	 * compute number of bits to split from sum
+	 */
+	dpsum = (pixelsum - (thisblock/2) - 1)/thisblock;
+	if (dpsum < 0) dpsum = 0.0;
+/*	psum = ((unsigned int) dpsum ) >> 1; */
+	psum = ((unsigned char) dpsum ) >> 1;
+	for (fs = 0; psum>0; fs++) psum >>= 1;
+
+	/*
+	 * write the codes
+	 * fsbits ID bits used to indicate split level
+	 */
+	if (fs >= fsmax) {
+	    /* Special high entropy case when FS >= fsmax
+	     * Just write pixel difference values directly, no Rice coding at all.
+	     */
+	    if (output_nbits(buffer, fsmax+1, fsbits) == EOF) {
+                ffpmsg("rice_encode: end of buffer");
+                free(diff);
+		return(-1);
+	    }
+	    for (j=0; j<thisblock; j++) {
+		if (output_nbits(buffer, diff[j], bbits) == EOF) {
+                    ffpmsg("rice_encode: end of buffer");
+                    free(diff);
+		    return(-1);
+		}
+	    }
+	} else if (fs == 0 && pixelsum == 0) {
+	    /*
+	     * special low entropy case when FS = 0 and pixelsum=0 (all
+	     * pixels in block are zero.)
+	     * Output a 0 and return
+	     */
+	    if (output_nbits(buffer, 0, fsbits) == EOF) {
+                ffpmsg("rice_encode: end of buffer");
+                free(diff);
+		return(-1);
+	    }
+	} else {
+	    /* normal case: not either very high or very low entropy */
+	    if (output_nbits(buffer, fs+1, fsbits) == EOF) {
+                ffpmsg("rice_encode: end of buffer");
+                free(diff);
+		return(-1);
+	    }
+	    fsmask = (1<<fs) - 1;
+	    /*
+	     * local copies of bit buffer to improve optimization
+	     */
+	    lbitbuffer = buffer->bitbuffer;
+	    lbits_to_go = buffer->bits_to_go;
+	    for (j=0; j<thisblock; j++) {
+		v = diff[j];
+		top = v >> fs;
+		/*
+		 * top is coded by top zeros + 1
+		 */
+		if (lbits_to_go >= top+1) {
+		    lbitbuffer <<= top+1;
+		    lbitbuffer |= 1;
+		    lbits_to_go -= top+1;
+		} else {
+		    lbitbuffer <<= lbits_to_go;
+		    putcbuf(lbitbuffer & 0xff,buffer);
+		    for (top -= lbits_to_go; top>=8; top -= 8) {
+			putcbuf(0, buffer);
+		    }
+		    lbitbuffer = 1;
+		    lbits_to_go = 7-top;
+		}
+		/*
+		 * bottom FS bits are written without coding
+		 * code is output_nbits, moved into this routine to reduce overheads
+		 * This code potentially breaks if FS>24, so I am limiting
+		 * FS to 24 by choice of FSMAX above.
+		 */
+		if (fs > 0) {
+		    lbitbuffer <<= fs;
+		    lbitbuffer |= v & fsmask;
+		    lbits_to_go -= fs;
+		    while (lbits_to_go <= 0) {
+			putcbuf((lbitbuffer>>(-lbits_to_go)) & 0xff,buffer);
+			lbits_to_go += 8;
+		    }
+		}
+	    }
+	    /* check if overflowed output buffer */
+	    if (buffer->current > buffer->end) {
+                 ffpmsg("rice_encode: end of buffer");
+                 free(diff);
+		 return(-1);
 	    }
 	    buffer->bitbuffer = lbitbuffer;
@@ -270,4 +743,11 @@
 int lbitbuffer;
 int lbits_to_go;
+    /* AND mask for the right-most n bits */
+    static unsigned int mask[33] = 
+         {0,
+	  0x1,       0x3,       0x7,       0xf,       0x1f,       0x3f,       0x7f,       0xff,
+	  0x1ff,     0x3ff,     0x7ff,     0xfff,     0x1fff,     0x3fff,     0x7fff,     0xffff,
+	  0x1ffff,   0x3ffff,   0x7ffff,   0xfffff,   0x1fffff,   0x3fffff,   0x7fffff,   0xffffff,
+	  0x1ffffff, 0x3ffffff, 0x7ffffff, 0xfffffff, 0x1fffffff, 0x3fffffff, 0x7fffffff, 0xffffffff};
 
     /*
@@ -282,11 +762,13 @@
 	 */
 	lbitbuffer <<= lbits_to_go;
-	lbitbuffer |= (bits>>(n-lbits_to_go)) & ((1<<lbits_to_go)-1);
-	if (putcbuf(lbitbuffer & 0xff,buffer) == EOF) return(EOF);
+/*	lbitbuffer |= (bits>>(n-lbits_to_go)) & ((1<<lbits_to_go)-1); */
+	lbitbuffer |= (bits>>(n-lbits_to_go)) & *(mask+lbits_to_go);
+	putcbuf(lbitbuffer & 0xff,buffer);
 	n -= lbits_to_go;
 	lbits_to_go = 8;
     }
     lbitbuffer <<= n;
-    lbitbuffer |= ( bits & ((1<<n)-1) );
+/*    lbitbuffer |= ( bits & ((1<<n)-1) ); */
+    lbitbuffer |= ( bits & *(mask+n) );
     lbits_to_go -= n;
     while (lbits_to_go <= 0) {
@@ -294,6 +776,5 @@
 	 * bitbuffer full, put out top 8 bits
 	 */
-	if (putcbuf((lbitbuffer>>(-lbits_to_go)) & 0xff,buffer) == EOF)
-	    return(EOF);
+	putcbuf((lbitbuffer>>(-lbits_to_go)) & 0xff,buffer);
 	lbits_to_go += 8;
     }
@@ -302,5 +783,4 @@
     return(0);
 }
-
 /*---------------------------------------------------------------------------*/
 /* Flush out the last bits */
@@ -309,6 +789,9 @@
 {
     if(buffer->bits_to_go < 8) {
-	if (putcbuf(buffer->bitbuffer<<buffer->bits_to_go,buffer) == EOF)
+	putcbuf(buffer->bitbuffer<<buffer->bits_to_go,buffer);
+	
+/*	if (putcbuf(buffer->bitbuffer<<buffer->bits_to_go,buffer) == EOF)
 	    return(EOF);
+*/
     }
     return(0);
@@ -334,4 +817,5 @@
 */
 
+/*---------------------------------------------------------------------------*/
 /* this routine used to be called 'rdecomp'  (WDP)  */
 
@@ -356,5 +840,6 @@
      */
     bsize = 4;
-/*    nblock = 32; */
+
+/*    nblock = 32; now an input parameter */
     /*
      * From bsize derive:
@@ -363,4 +848,6 @@
      * BBITS = bits/pixel for direct coding
      */
+
+/*
     switch (bsize) {
     case 1:
@@ -377,12 +864,15 @@
 	break;
     default:
-        fprintf (stderr, "rdecomp: bsize must be 1, 2, or 4 bytes");
+        ffpmsg("rdecomp: bsize must be 1, 2, or 4 bytes");
 	return 1;
     }
+*/
+
+    /* move out of switch block, to tweak performance */
+    fsbits = 5;
+    fsmax = 25;
+
     bbits = 1<<fsbits;
 
-    // XXX EAM : sensible code would put this in an init function,
-    // with a provided cleanup function.  this is just generating
-    // the number of 1 bits for any 8bit value
     if (nonzero_count == (int *) NULL) {
 	/*
@@ -394,5 +884,5 @@
 	nonzero_count = (int *) malloc(256*sizeof(int));
 	if (nonzero_count == (int *) NULL) {
-            fprintf (stderr, "rdecomp: insufficient memory");
+            ffpmsg("rdecomp: insufficient memory");
 	    return 1;
 	}
@@ -411,6 +901,4 @@
     /* first 4 bytes of input buffer contain the value of the first */
     /* 4 byte integer value, without any encoding */
-    
-    // XXX EAM : wow! this is just doing { lastpix = *(int *)c; } without a cast (encoded value is big-endian)
     
     lastpix = 0;
@@ -437,4 +925,5 @@
 	}
 	fs = (b >> nbits) - 1;
+
 	b &= (1<<nbits)-1;
 	/* loop over the next block */
@@ -494,4 +983,5 @@
 		diff = (nzero<<fs) | (b>>nbits);
 		b &= (1<<nbits)-1;
+
 		/* undo mapping and differencing */
 		if ((diff & 1) == 0) {
@@ -505,11 +995,378 @@
 	}
 	if (c > cend) {
-            fprintf (stderr, "decompression error: hit end of compressed byte stream");
+            ffpmsg("decompression error: hit end of compressed byte stream");
 	    return 1;
 	}
     }
     if (c < cend) {
-        fprintf (stderr, "decompression warning: unused bytes at end of compressed buffer");
+        ffpmsg("decompression warning: unused bytes at end of compressed buffer");
     }
     return 0;
 }
+/*---------------------------------------------------------------------------*/
+/* this routine used to be called 'rdecomp'  (WDP)  */
+
+int fits_rdecomp_short (unsigned char *c,		/* input buffer			*/
+	     int clen,			/* length of input		*/
+	     unsigned short array[],  	/* output array			*/
+	     int nx,			/* number of output pixels	*/
+	     int nblock)		/* coding block size		*/
+{
+int i, imax;
+int bsize, k;
+int nbits, nzero, fs;
+unsigned char *cend, bytevalue;
+unsigned int b, diff, lastpix;
+int fsmax, fsbits, bbits;
+static int *nonzero_count = (int *)NULL;
+
+   /*
+     * Original size of each pixel (bsize, bytes) and coding block
+     * size (nblock, pixels)
+     * Could make bsize a parameter to allow more efficient
+     * compression of short & byte images.
+     */
+
+    bsize = 2;
+    
+/*    nblock = 32; now an input parameter */
+    /*
+     * From bsize derive:
+     * FSBITS = # bits required to store FS
+     * FSMAX = maximum value for FS
+     * BBITS = bits/pixel for direct coding
+     */
+
+/*
+    switch (bsize) {
+    case 1:
+	fsbits = 3;
+	fsmax = 6;
+	break;
+    case 2:
+	fsbits = 4;
+	fsmax = 14;
+	break;
+    case 4:
+	fsbits = 5;
+	fsmax = 25;
+	break;
+    default:
+        ffpmsg("rdecomp: bsize must be 1, 2, or 4 bytes");
+	return 1;
+    }
+*/
+
+    /* move out of switch block, to tweak performance */
+    fsbits = 4;
+    fsmax = 14;
+
+    bbits = 1<<fsbits;
+
+    if (nonzero_count == (int *) NULL) {
+	/*
+	 * nonzero_count is lookup table giving number of bits
+	 * in 8-bit values not including leading zeros
+	 */
+
+        /*  NOTE!!!  This memory never gets freed  */
+	nonzero_count = (int *) malloc(256*sizeof(int));
+	if (nonzero_count == (int *) NULL) {
+            ffpmsg("rdecomp: insufficient memory");
+	    return 1;
+	}
+	nzero = 8;
+	k = 128;
+	for (i=255; i>=0; ) {
+	    for ( ; i>=k; i--) nonzero_count[i] = nzero;
+	    k = k/2;
+	    nzero--;
+	}
+    }
+    /*
+     * Decode in blocks of nblock pixels
+     */
+
+    /* first 2 bytes of input buffer contain the value of the first */
+    /* 2 byte integer value, without any encoding */
+    
+    lastpix = 0;
+    bytevalue = c[0];
+    lastpix = lastpix | (bytevalue<<8);
+    bytevalue = c[1];
+    lastpix = lastpix | bytevalue;
+
+    c += 2;  
+    cend = c + clen - 2;
+
+    b = *c++;		    /* bit buffer			*/
+    nbits = 8;		    /* number of bits remaining in b	*/
+    for (i = 0; i<nx; ) {
+	/* get the FS value from first fsbits */
+	nbits -= fsbits;
+	while (nbits < 0) {
+	    b = (b<<8) | (*c++);
+	    nbits += 8;
+	}
+	fs = (b >> nbits) - 1;
+
+	b &= (1<<nbits)-1;
+	/* loop over the next block */
+	imax = i + nblock;
+	if (imax > nx) imax = nx;
+	if (fs<0) {
+	    /* low-entropy case, all zero differences */
+	    for ( ; i<imax; i++) array[i] = lastpix;
+	} else if (fs==fsmax) {
+	    /* high-entropy case, directly coded pixel values */
+	    for ( ; i<imax; i++) {
+		k = bbits - nbits;
+		diff = b<<k;
+		for (k -= 8; k >= 0; k -= 8) {
+		    b = *c++;
+		    diff |= b<<k;
+		}
+		if (nbits>0) {
+		    b = *c++;
+		    diff |= b>>(-k);
+		    b &= (1<<nbits)-1;
+		} else {
+		    b = 0;
+		}
+   
+		/*
+		 * undo mapping and differencing
+		 * Note that some of these operations will overflow the
+		 * unsigned int arithmetic -- that's OK, it all works
+		 * out to give the right answers in the output file.
+		 */
+		if ((diff & 1) == 0) {
+		    diff = diff>>1;
+		} else {
+		    diff = ~(diff>>1);
+		}
+		array[i] = diff+lastpix;
+		lastpix = array[i];
+	    }
+	} else {
+	    /* normal case, Rice coding */
+	    for ( ; i<imax; i++) {
+		/* count number of leading zeros */
+		while (b == 0) {
+		    nbits += 8;
+		    b = *c++;
+		}
+		nzero = nbits - nonzero_count[b];
+		nbits -= nzero+1;
+		/* flip the leading one-bit */
+		b ^= 1<<nbits;
+		/* get the FS trailing bits */
+		nbits -= fs;
+		while (nbits < 0) {
+		    b = (b<<8) | (*c++);
+		    nbits += 8;
+		}
+		diff = (nzero<<fs) | (b>>nbits);
+		b &= (1<<nbits)-1;
+
+		/* undo mapping and differencing */
+		if ((diff & 1) == 0) {
+		    diff = diff>>1;
+		} else {
+		    diff = ~(diff>>1);
+		}
+		array[i] = diff+lastpix;
+		lastpix = array[i];
+	    }
+	}
+	if (c > cend) {
+            ffpmsg("decompression error: hit end of compressed byte stream");
+	    return 1;
+	}
+    }
+    if (c < cend) {
+        ffpmsg("decompression warning: unused bytes at end of compressed buffer");
+    }
+    return 0;
+}
+/*---------------------------------------------------------------------------*/
+/* this routine used to be called 'rdecomp'  (WDP)  */
+
+int fits_rdecomp_byte (unsigned char *c,		/* input buffer			*/
+	     int clen,			/* length of input		*/
+	     unsigned char array[],  	/* output array			*/
+	     int nx,			/* number of output pixels	*/
+	     int nblock)		/* coding block size		*/
+{
+int i, imax;
+int bsize, k;
+int nbits, nzero, fs;
+unsigned char *cend;
+unsigned int b, diff, lastpix;
+int fsmax, fsbits, bbits;
+static int *nonzero_count = (int *)NULL;
+
+   /*
+     * Original size of each pixel (bsize, bytes) and coding block
+     * size (nblock, pixels)
+     * Could make bsize a parameter to allow more efficient
+     * compression of short & byte images.
+     */
+
+    bsize = 1;
+    
+/*    nblock = 32; now an input parameter */
+    /*
+     * From bsize derive:
+     * FSBITS = # bits required to store FS
+     * FSMAX = maximum value for FS
+     * BBITS = bits/pixel for direct coding
+     */
+
+/*
+    switch (bsize) {
+    case 1:
+	fsbits = 3;
+	fsmax = 6;
+	break;
+    case 2:
+	fsbits = 4;
+	fsmax = 14;
+	break;
+    case 4:
+	fsbits = 5;
+	fsmax = 25;
+	break;
+    default:
+        ffpmsg("rdecomp: bsize must be 1, 2, or 4 bytes");
+	return 1;
+    }
+*/
+
+    /* move out of switch block, to tweak performance */
+    fsbits = 3;
+    fsmax = 6;
+
+    bbits = 1<<fsbits;
+
+    if (nonzero_count == (int *) NULL) {
+	/*
+	 * nonzero_count is lookup table giving number of bits
+	 * in 8-bit values not including leading zeros
+	 */
+
+        /*  NOTE!!!  This memory never gets freed  */
+	nonzero_count = (int *) malloc(256*sizeof(int));
+	if (nonzero_count == (int *) NULL) {
+            ffpmsg("rdecomp: insufficient memory");
+	    return 1;
+	}
+	nzero = 8;
+	k = 128;
+	for (i=255; i>=0; ) {
+	    for ( ; i>=k; i--) nonzero_count[i] = nzero;
+	    k = k/2;
+	    nzero--;
+	}
+    }
+    /*
+     * Decode in blocks of nblock pixels
+     */
+
+    /* first byte of input buffer contain the value of the first */
+    /* byte integer value, without any encoding */
+    
+    lastpix = c[0];
+    c += 1;  
+    cend = c + clen - 1;
+
+    b = *c++;		    /* bit buffer			*/
+    nbits = 8;		    /* number of bits remaining in b	*/
+    for (i = 0; i<nx; ) {
+	/* get the FS value from first fsbits */
+	nbits -= fsbits;
+	while (nbits < 0) {
+	    b = (b<<8) | (*c++);
+	    nbits += 8;
+	}
+	fs = (b >> nbits) - 1;
+
+	b &= (1<<nbits)-1;
+	/* loop over the next block */
+	imax = i + nblock;
+	if (imax > nx) imax = nx;
+	if (fs<0) {
+	    /* low-entropy case, all zero differences */
+	    for ( ; i<imax; i++) array[i] = lastpix;
+	} else if (fs==fsmax) {
+	    /* high-entropy case, directly coded pixel values */
+	    for ( ; i<imax; i++) {
+		k = bbits - nbits;
+		diff = b<<k;
+		for (k -= 8; k >= 0; k -= 8) {
+		    b = *c++;
+		    diff |= b<<k;
+		}
+		if (nbits>0) {
+		    b = *c++;
+		    diff |= b>>(-k);
+		    b &= (1<<nbits)-1;
+		} else {
+		    b = 0;
+		}
+   
+		/*
+		 * undo mapping and differencing
+		 * Note that some of these operations will overflow the
+		 * unsigned int arithmetic -- that's OK, it all works
+		 * out to give the right answers in the output file.
+		 */
+		if ((diff & 1) == 0) {
+		    diff = diff>>1;
+		} else {
+		    diff = ~(diff>>1);
+		}
+		array[i] = diff+lastpix;
+		lastpix = array[i];
+	    }
+	} else {
+	    /* normal case, Rice coding */
+	    for ( ; i<imax; i++) {
+		/* count number of leading zeros */
+		while (b == 0) {
+		    nbits += 8;
+		    b = *c++;
+		}
+		nzero = nbits - nonzero_count[b];
+		nbits -= nzero+1;
+		/* flip the leading one-bit */
+		b ^= 1<<nbits;
+		/* get the FS trailing bits */
+		nbits -= fs;
+		while (nbits < 0) {
+		    b = (b<<8) | (*c++);
+		    nbits += 8;
+		}
+		diff = (nzero<<fs) | (b>>nbits);
+		b &= (1<<nbits)-1;
+
+		/* undo mapping and differencing */
+		if ((diff & 1) == 0) {
+		    diff = diff>>1;
+		} else {
+		    diff = ~(diff>>1);
+		}
+		array[i] = diff+lastpix;
+		lastpix = array[i];
+	    }
+	}
+	if (c > cend) {
+            ffpmsg("decompression error: hit end of compressed byte stream");
+	    return 1;
+	}
+    }
+    if (c < cend) {
+        ffpmsg("decompression warning: unused bytes at end of compressed buffer");
+    }
+    return 0;
+}
Index: /trunk/Ohana/src/libfits/extern/ricecomp.h
===================================================================
--- /trunk/Ohana/src/libfits/extern/ricecomp.h	(revision 18268)
+++ /trunk/Ohana/src/libfits/extern/ricecomp.h	(revision 18269)
@@ -1,107 +1,12 @@
-/* @(#) buffer.h 1.1 98/07/21 12:34:27 */
-/* buffer.h: structure for compression to buffer rather than to a file, including
- *   bit I/O buffer
- *
- * R. White, 19 June 1998
- */
 
 
-typedef unsigned char Buffer_t;
+// 2008.06.22 EAM : This code was taken from CFITSIO and included in Ohana for rice
+// decompression.  In order to include this .c file (or future upgrades), we include our own
+// version of ricecomp.h (not the one available in the CFITSIO tree) and include it instead of
+// fitsio2.h.  This file defines our own version of ffpmsg.  
 
-typedef struct {
-	int bitbuffer;		/* bit buffer					*/
-	int bits_to_go;		/* bits to go in buffer			*/
-	Buffer_t *start;	/* start of buffer				*/
-	Buffer_t *current;	/* current position in buffer	*/
-	Buffer_t *end;		/* end of buffer				*/
-} Buffer;
-
-#define buffree(mf)		(free(mf->start), free(mf))
-#define bufused(mf)		(mf->current - mf->start)
-#define bufreset(mf)	(mf->current = mf->start)
-
-/*
- * getcbuf, putcbuf macros for character IO to buffer
- * putcbuf returns EOF on end of buffer, else returns 0
- */
-#define getcbuf(mf) ((mf->current >= mf->end) ? EOF : *(mf->current)++)
-#define putcbuf(c,mf) \
-	((mf->current >= mf->end) ? \
-	EOF :\
-	((*(mf->current)++ = c), 0))
-
-/*
- * bufalloc sets up buffer of length n
- */
-
-/*  not needed by CFITSIO
-
-static Buffer *bufalloc(int n)
-{
-Buffer *mf;
-
-	mf = (Buffer *) malloc(sizeof(Buffer));
-	if (mf == (Buffer *)NULL) return((Buffer *)NULL);
-
-	mf->start = (Buffer_t *) malloc(n*sizeof(Buffer_t));
-	if (mf->start == (Buffer_t *)NULL) {
-		free(mf);
-		return((Buffer *)NULL);
-	}
-	mf->bits_to_go = 8;
-	mf->end = mf->start + n;
-	mf->current = mf->start;
-	return(mf);
+int ffpmsg (char *line) {
+    
+    fprintf (stderr, "%s\n", line);
+    return 1;
 }
-*/
-
-/*
- * bufrealloc extends buffer (or truncates it) by
- * reallocating memory
- */
-
-/* not needed by CFITSIO 
-static int bufrealloc(Buffer *mf, int n)
-{
-int len;
-
-	len = mf->current - mf->start;
-
-	* silently throw away data if buffer is already longer than n *
-	if (len>n) len = n;
-	if (len<0) len = 0;
-
-	mf->start = (Buffer_t *) realloc(mf->start, n*sizeof(Buffer_t));
-	if (mf->start == (Buffer_t *)NULL) return(0);
-
-	mf->end = mf->start + n;
-	mf->current = mf->start + len;
-	return(n);
-}
-*/
-
-/*
- * bufdump dumps contents of buffer to outfile and resets
- * it to be empty.  Returns number of bytes written.
- *
- * Note we don't write out the bit buffer -- you must call
- * done_outputing_bits() first to ensure that the bit buffer
- * is written out.  I do it this way to allow incremental
- * buffer dumps while bit IO is still going on.
- */
-
-/*  not needed by CFITSIO
-
-static int bufdump(FILE *outfile, Buffer *buffer)
-{
-int ndump;
-
-	ndump = bufused(buffer);
-	if (fwrite(buffer->start, 1, ndump, outfile) != ndump) {
-		fprintf(stderr, "bufdump: error in write\n");
-		exit(1);
-    }
-	bufreset(buffer);
-	return(ndump);
-}
-*/