« back to all changes in this revision
Viewing changes to ricecomp.c
- browse files at revision 1
- compare with another revision
- download diff
- download tarball
- view history from revision 1
- Committer: Bazaar Package Importer
- Author(s): Gopal Narayanan
- Date: 2002-02-26 11:27:29 UTC
- Revision ID: james.westby@ubuntu.com-20020226112729-3q2o993rhh81ipp4
Tags: upstream-2.401
Import upstream version 2.401
- files added:
added
removed
ricecomp.c
1
/*
2
The following code was written by Richard White at STScI and made
3
available for use in CFITSIO in July 1999. These routines were
4
originally contained in 2 source files: rcomp.c and rdecomp.c,
5
and the 'include' file now called ricecomp.h was originally called buffer.h.
6
*/
7
8
/*----------------------------------------------------------*/
9
/* */
10
/* START OF SOURCE FILE ORIGINALLY CALLED rcomp.c */
11
/* */
12
/*----------------------------------------------------------*/
13
/* @(#) rcomp.c 1.5 99/03/01 12:40:27 */
14
/* rcomp.c Compress image line using
15
* (1) Difference of adjacent pixels
16
* (2) Rice algorithm coding
17
*
18
* Returns number of bytes written to code buffer or
19
* -1 on failure
20
*/
21
22
#include <stdio.h>
23
#include <stdlib.h>
24
#include <string.h>
25
#include "ricecomp.h" /* originally included in rcomp.c file (WDP) */
26
#include "fitsio2.h"
27
28
29
static void start_outputing_bits(Buffer *buffer);
30
static int done_outputing_bits(Buffer *buffer);
31
static int output_nbits(Buffer *buffer, int bits, int n);
32
33
/* this routine used to be called 'rcomp' (WDP) */
34
35
int fits_rcomp(int a[], /* input array */
36
int nx, /* number of input pixels */
37
unsigned char *c, /* output buffer */
38
int clen, /* max length of output */
39
int nblock) /* coding block size */
40
{
41
Buffer bufmem, *buffer = &bufmem;
42
int bsize, i, j, thisblock;
43
int lastpix, nextpix, pdiff;
44
int v, fs, fsmask, top, fsmax, fsbits, bbits;
45
int lbitbuffer, lbits_to_go;
46
unsigned int psum;
47
double pixelsum, dpsum;
48
unsigned int *diff;
49
50
/*
51
* Original size of each pixel (bsize, bytes) and coding block
52
* size (nblock, pixels)
53
* Could make bsize a parameter to allow more efficient
54
* compression of short & byte images.
55
*/
56
bsize = 4;
57
/* nblock = 32; */
58
/*
59
* From bsize derive:
60
* FSBITS = # bits required to store FS
61
* FSMAX = maximum value for FS
62
* BBITS = bits/pixel for direct coding
63
*/
64
switch (bsize) {
65
case 1:
66
fsbits = 3;
67
fsmax = 6;
68
break;
69
case 2:
70
fsbits = 4;
71
fsmax = 14;
72
break;
73
case 4:
74
fsbits = 5;
75
fsmax = 25;
76
break;
77
default:
78
ffpmsg("rdecomp: bsize must be 1, 2, or 4 bytes");
79
return(-1);
80
}
81
bbits = 1<<fsbits;
82
83
/*
84
* Set up buffer pointers
85
*/
86
buffer->start = c;
87
buffer->current = c;
88
buffer->end = c+clen;
89
buffer->bits_to_go = 8;
90
/*
91
* array for differences mapped to non-negative values
92
*/
93
diff = (unsigned int *) malloc(nblock*sizeof(unsigned int));
94
if (diff == (unsigned int *) NULL) {
95
ffpmsg("fits_rcomp: insufficient memory");
96
return(-1);
97
}
98
/*
99
* Code in blocks of nblock pixels
100
*/
101
start_outputing_bits(buffer);
102
103
/* write out first int value to the first 4 bytes of the buffer */
104
if (output_nbits(buffer, a[0], 32) == EOF) {
105
ffpmsg("rice_encode: end of buffer");
106
free(diff);
107
return(-1);
108
}
109
110
lastpix = a[0]; /* the first difference will always be zero */
111
112
thisblock = nblock;
113
for (i=0; i<nx; i += nblock) {
114
/* last block may be shorter */
115
if (nx-i < nblock) thisblock = nx-i;
116
/*
117
* Compute differences of adjacent pixels and map them to unsigned values.
118
* Note that this may overflow the integer variables -- that's
119
* OK, because we can recover when decompressing. If we were
120
* compressing shorts or bytes, would want to do this arithmetic
121
* with short/byte working variables (though diff will still be
122
* passed as an int.)
123
*
124
* compute sum of mapped pixel values at same time
125
* use double precision for sum to allow 32-bit integer inputs
126
*/
127
pixelsum = 0.0;
128
for (j=0; j<thisblock; j++) {
129
nextpix = a[i+j];
130
pdiff = nextpix - lastpix;
131
diff[j] = (unsigned int) ((pdiff<0) ? ~(pdiff<<1) : (pdiff<<1));
132
pixelsum += diff[j];
133
lastpix = nextpix;
134
}
135
/*
136
* compute number of bits to split from sum
137
*/
138
dpsum = (pixelsum - (thisblock/2) - 1)/thisblock;
139
if (dpsum < 0) dpsum = 0.0;
140
psum = ((unsigned int) dpsum ) >> 1;
141
for (fs = 0; psum>0; fs++) psum >>= 1;
142
/*
143
* write the codes
144
* fsbits ID bits used to indicate split level
145
*/
146
if (fs >= fsmax) {
147
/* Special high entropy case when FS >= fsmax
148
* Just write pixel difference values directly, no Rice coding at all.
149
*/
150
if (output_nbits(buffer, fsmax+1, fsbits) == EOF) {
151
ffpmsg("rice_encode: end of buffer");
152
free(diff);
153
return(-1);
154
}
155
for (j=0; j<thisblock; j++) {
156
if (output_nbits(buffer, diff[j], bbits) == EOF) {
157
ffpmsg("rice_encode: end of buffer");
158
free(diff);
159
return(-1);
160
}
161
}
162
} else if (fs == 0 && pixelsum == 0) {
163
/*
164
* special low entropy case when FS = 0 and pixelsum=0 (all
165
* pixels in block are zero.)
166
* Output a 0 and return
167
*/
168
if (output_nbits(buffer, 0, fsbits) == EOF) {
169
ffpmsg("rice_encode: end of buffer");
170
free(diff);
171
return(-1);
172
}
173
} else {
174
/* normal case: not either very high or very low entropy */
175
if (output_nbits(buffer, fs+1, fsbits) == EOF) {
176
ffpmsg("rice_encode: end of buffer");
177
free(diff);
178
return(-1);
179
}
180
fsmask = (1<<fs) - 1;
181
/*
182
* local copies of bit buffer to improve optimization
183
*/
184
lbitbuffer = buffer->bitbuffer;
185
lbits_to_go = buffer->bits_to_go;
186
for (j=0; j<thisblock; j++) {
187
v = diff[j];
188
top = v >> fs;
189
/*
190
* top is coded by top zeros + 1
191
*/
192
if (lbits_to_go >= top+1) {
193
lbitbuffer <<= top+1;
194
lbitbuffer |= 1;
195
lbits_to_go -= top+1;
196
} else {
197
lbitbuffer <<= lbits_to_go;
198
if (putcbuf(lbitbuffer & 0xff,buffer) == EOF) {
199
ffpmsg("rice_encode: end of buffer");
200
free(diff);
201
return(-1);
202
}
203
for (top -= lbits_to_go; top>=8; top -= 8) {
204
if (putcbuf(0, buffer) == EOF) {
205
ffpmsg("rice_encode: end of buffer");
206
free(diff);
207
return(-1);
208
}
209
}
210
lbitbuffer = 1;
211
lbits_to_go = 7-top;
212
}
213
/*
214
* bottom FS bits are written without coding
215
* code is output_nbits, moved into this routine to reduce overheads
216
* This code potentially breaks if FS>24, so I am limiting
217
* FS to 24 by choice of FSMAX above.
218
*/
219
if (fs > 0) {
220
lbitbuffer <<= fs;
221
lbitbuffer |= v & fsmask;
222
lbits_to_go -= fs;
223
while (lbits_to_go <= 0) {
224
if (putcbuf((lbitbuffer>>(-lbits_to_go)) & 0xff,buffer)==EOF) {
225
ffpmsg("rice_encode: end of buffer");
226
free(diff);
227
return(-1);
228
}
229
lbits_to_go += 8;
230
}
231
}
232
}
233
buffer->bitbuffer = lbitbuffer;
234
buffer->bits_to_go = lbits_to_go;
235
}
236
}
237
done_outputing_bits(buffer);
238
free(diff);
239
/*
240
* return number of bytes used
241
*/
242
return(buffer->current - buffer->start);
243
}
244
/*---------------------------------------------------------------------------*/
245
/* bit_output.c
246
*
247
* Bit output routines
248
* Procedures return zero on success, EOF on end-of-buffer
249
*
250
* Programmer: R. White Date: 20 July 1998
251
*/
252
253
/* Initialize for bit output */
254
255
static void start_outputing_bits(Buffer *buffer)
256
{
257
/*
258
* Buffer is empty to start with
259
*/
260
buffer->bitbuffer = 0;
261
buffer->bits_to_go = 8;
262
}
263
264
/*---------------------------------------------------------------------------*/
265
/* Output N bits (N must be <= 32) */
266
267
static int output_nbits(Buffer *buffer, int bits, int n)
268
{
269
/* local copies */
270
int lbitbuffer;
271
int lbits_to_go;
272
273
/*
274
* insert bits at end of bitbuffer
275
*/
276
lbitbuffer = buffer->bitbuffer;
277
lbits_to_go = buffer->bits_to_go;
278
if (lbits_to_go+n > 32) {
279
/*
280
* special case for large n: put out the top lbits_to_go bits first
281
* note that 0 < lbits_to_go <= 8
282
*/
283
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);
286
n -= lbits_to_go;
287
lbits_to_go = 8;
288
}
289
lbitbuffer <<= n;
290
lbitbuffer |= ( bits & ((1<<n)-1) );
291
lbits_to_go -= n;
292
while (lbits_to_go <= 0) {
293
/*
294
* bitbuffer full, put out top 8 bits
295
*/
296
if (putcbuf((lbitbuffer>>(-lbits_to_go)) & 0xff,buffer) == EOF)
297
return(EOF);
298
lbits_to_go += 8;
299
}
300
buffer->bitbuffer = lbitbuffer;
301
buffer->bits_to_go = lbits_to_go;
302
return(0);
303
}
304
305
/*---------------------------------------------------------------------------*/
306
/* Flush out the last bits */
307
308
static int done_outputing_bits(Buffer *buffer)
309
{
310
if(buffer->bits_to_go < 8) {
311
if (putcbuf(buffer->bitbuffer<<buffer->bits_to_go,buffer) == EOF)
312
return(EOF);
313
}
314
return(0);
315
}
316
/*---------------------------------------------------------------------------*/
317
/*----------------------------------------------------------*/
318
/* */
319
/* START OF SOURCE FILE ORIGINALLY CALLED rdecomp.c */
320
/* */
321
/*----------------------------------------------------------*/
322
323
/* @(#) rdecomp.c 1.4 99/03/01 12:38:41 */
324
/* rdecomp.c Decompress image line using
325
* (1) Difference of adjacent pixels
326
* (2) Rice algorithm coding
327
*
328
* Returns 0 on success or 1 on failure
329
*/
330
331
/* moved these 'includes' to the beginning of the file (WDP)
332
#include <stdio.h>
333
#include <stdlib.h>
334
*/
335
336
/* this routine used to be called 'rdecomp' (WDP) */
337
338
int fits_rdecomp (unsigned char *c, /* input buffer */
339
int clen, /* length of input */
340
unsigned int array[], /* output array */
341
int nx, /* number of output pixels */
342
int nblock) /* coding block size */
343
{
344
int bsize, i, k, imax;
345
int nbits, nzero, fs;
346
unsigned char *cend, bytevalue;
347
unsigned int b, diff, lastpix;
348
int fsmax, fsbits, bbits;
349
static int *nonzero_count = (int *)NULL;
350
351
/*
352
* Original size of each pixel (bsize, bytes) and coding block
353
* size (nblock, pixels)
354
* Could make bsize a parameter to allow more efficient
355
* compression of short & byte images.
356
*/
357
bsize = 4;
358
/* nblock = 32; */
359
/*
360
* From bsize derive:
361
* FSBITS = # bits required to store FS
362
* FSMAX = maximum value for FS
363
* BBITS = bits/pixel for direct coding
364
*/
365
switch (bsize) {
366
case 1:
367
fsbits = 3;
368
fsmax = 6;
369
break;
370
case 2:
371
fsbits = 4;
372
fsmax = 14;
373
break;
374
case 4:
375
fsbits = 5;
376
fsmax = 25;
377
break;
378
default:
379
ffpmsg("rdecomp: bsize must be 1, 2, or 4 bytes");
380
return 1;
381
}
382
bbits = 1<<fsbits;
383
384
if (nonzero_count == (int *) NULL) {
385
/*
386
* nonzero_count is lookup table giving number of bits
387
* in 8-bit values not including leading zeros
388
*/
389
390
/* NOTE!!! This memory never gets freed */
391
nonzero_count = (int *) malloc(256*sizeof(int));
392
if (nonzero_count == (int *) NULL) {
393
ffpmsg("rdecomp: insufficient memory");
394
return 1;
395
}
396
nzero = 8;
397
k = 128;
398
for (i=255; i>=0; ) {
399
for ( ; i>=k; i--) nonzero_count[i] = nzero;
400
k = k/2;
401
nzero--;
402
}
403
}
404
/*
405
* Decode in blocks of nblock pixels
406
*/
407
408
/* first 4 bytes of input buffer contain the value of the first */
409
/* 4 byte integer value, without any encoding */
410
411
lastpix = 0;
412
bytevalue = c[0];
413
lastpix = lastpix | (bytevalue<<24);
414
bytevalue = c[1];
415
lastpix = lastpix | (bytevalue<<16);
416
bytevalue = c[2];
417
lastpix = lastpix | (bytevalue<<8);
418
bytevalue = c[3];
419
lastpix = lastpix | bytevalue;
420
421
c += 4;
422
cend = c + clen - 4;
423
424
b = *c++; /* bit buffer */
425
nbits = 8; /* number of bits remaining in b */
426
for (i = 0; i<nx; ) {
427
/* get the FS value from first fsbits */
428
nbits -= fsbits;
429
while (nbits < 0) {
430
b = (b<<8) | (*c++);
431
nbits += 8;
432
}
433
fs = (b >> nbits) - 1;
434
b &= (1<<nbits)-1;
435
/* loop over the next block */
436
imax = i + nblock;
437
if (imax > nx) imax = nx;
438
if (fs<0) {
439
/* low-entropy case, all zero differences */
440
for ( ; i<imax; i++) array[i] = lastpix;
441
} else if (fs==fsmax) {
442
/* high-entropy case, directly coded pixel values */
443
for ( ; i<imax; i++) {
444
k = bbits - nbits;
445
diff = b<<k;
446
for (k -= 8; k >= 0; k -= 8) {
447
b = *c++;
448
diff |= b<<k;
449
}
450
if (nbits>0) {
451
b = *c++;
452
diff |= b>>(-k);
453
b &= (1<<nbits)-1;
454
} else {
455
b = 0;
456
}
457
/*
458
* undo mapping and differencing
459
* Note that some of these operations will overflow the
460
* unsigned int arithmetic -- that's OK, it all works
461
* out to give the right answers in the output file.
462
*/
463
if ((diff & 1) == 0) {
464
diff = diff>>1;
465
} else {
466
diff = ~(diff>>1);
467
}
468
array[i] = diff+lastpix;
469
lastpix = array[i];
470
}
471
} else {
472
/* normal case, Rice coding */
473
for ( ; i<imax; i++) {
474
/* count number of leading zeros */
475
while (b == 0) {
476
nbits += 8;
477
b = *c++;
478
}
479
nzero = nbits - nonzero_count[b];
480
nbits -= nzero+1;
481
/* flip the leading one-bit */
482
b ^= 1<<nbits;
483
/* get the FS trailing bits */
484
nbits -= fs;
485
while (nbits < 0) {
486
b = (b<<8) | (*c++);
487
nbits += 8;
488
}
489
diff = (nzero<<fs) | (b>>nbits);
490
b &= (1<<nbits)-1;
491
/* undo mapping and differencing */
492
if ((diff & 1) == 0) {
493
diff = diff>>1;
494
} else {
495
diff = ~(diff>>1);
496
}
497
array[i] = diff+lastpix;
498
lastpix = array[i];
499
}
500
}
501
if (c > cend) {
502
ffpmsg("decompression error: hit end of compressed byte stream");
503
return 1;
504
}
505
}
506
if (c < cend) {
507
ffpmsg("decompression warning: unused bytes at end of compressed buffer");
508
}
509
return 0;
510
}
Loggerhead is a web-based interface for Breezy
Version: 2.0.1