« back to all changes in this revision
Viewing changes to src/fits_hcompress.c
- browse files at revision 17
- compare with another revision
- download diff
- download tarball
- view history from revision 17
- Committer: Package Import Robot
- Author(s): Aurelien Jarno
- Date: 2013-12-07 16:18:48 UTC
- mfrom: (1.1.11)
- Revision ID: package-import@ubuntu.com-20131207161848-mcw0saz0iprjhbju
Tags: 1:3.2-1
* New upstream version.
* Bump Standards-Version to 3.9.5 (no changes).
* Remove build-depends on zlib1g-dev and remove patches/01-zlib.diff.
* Add build-depends on libcfitsio3-dev and add
patches/01-system-cfitsio.diff.
* Update debian/copyright.
* Install upstream changelog now that it is provided in the upstream
tarball.
* Don't compress the binary packages with xz, it's no the dpkg's default.
* Bump Standards-Version to 3.9.5 (no changes).
* Remove build-depends on zlib1g-dev and remove patches/01-zlib.diff.
* Add build-depends on libcfitsio3-dev and add
patches/01-system-cfitsio.diff.
* Update debian/copyright.
* Install upstream changelog now that it is provided in the upstream
tarball.
* Don't compress the binary packages with xz, it's no the dpkg's default.
- files added:
- .pc/01-system-cfitsio.diff
- cextern
- cextern/cfitsio
- files removed:
- .pc/01-zlib.diff
- .pc/01-zlib.diff/src
- files modified:
- .pc/applied-patches
added
removed
src/fits_hcompress.c
1
/* $Id$
2
*/
3
4
/*****************************************************************************/
5
/* */
6
/* This file, fits_hcompress.c, contains the code required to compress */
7
/* data using the HCOMPRESS_1 compression format. */
8
/* */
9
/* Copyright (C) 2004 Association of Universities for Research in Astronomy */
10
/* (AURA) */
11
/* */
12
/* Redistribution and use in source and binary forms, with or without */
13
/* modification, are permitted provided that the following conditions are */
14
/* met: */
15
/* */
16
/* 1. Redistributions of source code must retain the above copyright */
17
/* notice, this list of conditions and the following disclaimer. */
18
/* */
19
/* 2. Redistributions in binary form must reproduce the above */
20
/* copyright notice, this list of conditions and the following */
21
/* disclaimer in the documentation and/or other materials provided */
22
/* with the distribution. */
23
/* */
24
/* 3. The name of AURA and its representatives may not be used to */
25
/* endorse or promote products derived from this software without */
26
/* specific prior written permission. */
27
/* */
28
/* THIS SOFTWARE IS PROVIDED BY AURA ``AS IS'' AND ANY EXPRESS OR IMPLIED */
29
/* WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF */
30
/* MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE */
31
/* DISCLAIMED. IN NO EVENT SHALL AURA BE LIABLE FOR ANY DIRECT, INDIRECT, */
32
/* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, */
33
/* BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS */
34
/* OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND */
35
/* ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR */
36
/* TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE */
37
/* USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH */
38
/* DAMAGE. */
39
/* */
40
/* This file was copied intact from the FITSIO software that was written by */
41
/* William Pence at the High Energy Astrophysic Science Archive Research */
42
/* Center (HEASARC) at the NASA Goddard Space Flight Center. That software */
43
/* contained the following copyright and warranty notices: */
44
/* */
45
/* Copyright (Unpublished--all rights reserved under the copyright laws of */
46
/* the United States), U.S. Government as represented by the Administrator */
47
/* of the National Aeronautics and Space Administration. No copyright is */
48
/* claimed in the United States under Title 17, U.S. Code. */
49
/* */
50
/* Permission to freely use, copy, modify, and distribute this software */
51
/* and its documentation without fee is hereby granted, provided that this */
52
/* copyright notice and disclaimer of warranty appears in all copies. */
53
/* */
54
/* DISCLAIMER: */
55
/* */
56
/* THE SOFTWARE IS PROVIDED 'AS IS' WITHOUT ANY WARRANTY OF ANY KIND, */
57
/* EITHER EXPRESSED, IMPLIED, OR STATUTORY, INCLUDING, BUT NOT LIMITED TO, */
58
/* ANY WARRANTY THAT THE SOFTWARE WILL CONFORM TO SPECIFICATIONS, ANY */
59
/* IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR */
60
/* PURPOSE, AND FREEDOM FROM INFRINGEMENT, AND ANY WARRANTY THAT THE */
61
/* DOCUMENTATION WILL CONFORM TO THE SOFTWARE, OR ANY WARRANTY THAT THE */
62
/* SOFTWARE WILL BE ERROR FREE. IN NO EVENT SHALL NASA BE LIABLE FOR ANY */
63
/* DAMAGES, INCLUDING, BUT NOT LIMITED TO, DIRECT, INDIRECT, SPECIAL OR */
64
/* CONSEQUENTIAL DAMAGES, ARISING OUT OF, RESULTING FROM, OR IN ANY WAY */
65
/* CONNECTED WITH THIS SOFTWARE, WHETHER OR NOT BASED UPON WARRANTY, */
66
/* CONTRACT, TORT , OR OTHERWISE, WHETHER OR NOT INJURY WAS SUSTAINED BY */
67
/* PERSONS OR PROPERTY OR OTHERWISE, AND WHETHER OR NOT LOSS WAS SUSTAINED */
68
/* FROM, OR AROSE OUT OF THE RESULTS OF, OR USE OF, THE SOFTWARE OR */
69
/* SERVICES PROVIDED HEREUNDER." */
70
/* */
71
/*****************************************************************************/
72
73
/* #########################################################################
74
These routines to apply the H-compress compression algorithm to a 2-D Fits
75
image were written by R. White at the STScI and were obtained from the STScI at
76
http://www.stsci.edu/software/hcompress.html
77
78
This source file is a concatination of the following sources files in the
79
original distribution
80
htrans.c
81
digitize.c
82
encode.c
83
qwrite.c
84
doencode.c
85
bit_output.c
86
qtree_encode.c
87
88
The following modifications have been made to the original code:
89
90
- commented out redundant "include" statements
91
- added the noutchar global variable
92
- changed all the 'extern' declarations to 'static', since all the routines are in
93
the same source file
94
- changed the first parameter in encode (and in lower level routines from a file stream
95
to a char array
96
- modifid the encode routine to return the size of the compressed array of bytes
97
- changed calls to printf and perror to call the CFITSIO ffpmsg routine
98
- modified the mywrite routine, and lower level byte writing routines, to copy
99
the output bytes to a char array, instead of writing them to a file stream
100
- replace "exit" statements with "return" statements
101
- changed the function declarations to the more modern ANSI C style
102
103
############################################################################ */
104
105
#include <stdio.h>
106
#include <string.h>
107
#include <math.h>
108
#include <stdlib.h>
109
#include "fitsio.h"
110
111
static long noutchar;
112
static long noutmax;
113
114
static int htrans(int a[],int nx,int ny);
115
static void digitize(int a[], int nx, int ny, int scale);
116
static int encode(char *outfile, long *nlen, int a[], int nx, int ny, int scale);
117
static void shuffle(int a[], int n, int n2, int tmp[]);
118
119
static int htrans64(LONGLONG a[],int nx,int ny);
120
static void digitize64(LONGLONG a[], int nx, int ny, int scale);
121
static int encode64(char *outfile, long *nlen, LONGLONG a[], int nx, int ny, int scale);
122
static void shuffle64(LONGLONG a[], int n, int n2, LONGLONG tmp[]);
123
124
static void writeint(char *outfile, int a);
125
static void writelonglong(char *outfile, LONGLONG a);
126
static int doencode(char *outfile, int a[], int nx, int ny, unsigned char nbitplanes[3]);
127
static int doencode64(char *outfile, LONGLONG a[], int nx, int ny, unsigned char nbitplanes[3]);
128
static int qwrite(char *file, char buffer[], int n);
129
130
static int qtree_encode(char *outfile, int a[], int n, int nqx, int nqy, int nbitplanes);
131
static int qtree_encode64(char *outfile, LONGLONG a[], int n, int nqx, int nqy, int nbitplanes);
132
static void start_outputing_bits(void);
133
static void done_outputing_bits(char *outfile);
134
static void output_nbits(char *outfile, int bits, int n);
135
136
static void qtree_onebit(int a[], int n, int nx, int ny, unsigned char b[], int bit);
137
static void qtree_onebit64(LONGLONG a[], int n, int nx, int ny, unsigned char b[], int bit);
138
static void qtree_reduce(unsigned char a[], int n, int nx, int ny, unsigned char b[]);
139
static int bufcopy(unsigned char a[], int n, unsigned char buffer[], int *b, int bmax);
140
static void write_bdirect(char *outfile, int a[], int n,int nqx, int nqy, unsigned char scratch[], int bit);
141
static void write_bdirect64(char *outfile, LONGLONG a[], int n,int nqx, int nqy, unsigned char scratch[], int bit);
142
143
/* #define output_nybble(outfile,c) output_nbits(outfile,c,4) */
144
static void output_nybble(char *outfile, int bits);
145
static void output_nnybble(char *outfile, int n, unsigned char array[]);
146
147
#define output_huffman(outfile,c) output_nbits(outfile,code[c],ncode[c])
148
149
/* ---------------------------------------------------------------------- */
150
int _pyfits_fits_hcompress(int *a, int ny, int nx, int scale, char *output,
151
long *nbytes, int *status)
152
{
153
/*
154
compress the input image using the H-compress algorithm
155
156
a - input image array
157
nx - size of X axis of image
158
ny - size of Y axis of image
159
scale - quantization scale factor. Larger values results in more (lossy) compression
160
scale = 0 does lossless compression
161
output - pre-allocated array to hold the output compressed stream of bytes
162
nbyts - input value = size of the output buffer;
163
returned value = size of the compressed byte stream, in bytes
164
165
NOTE: the nx and ny dimensions as defined within this code are reversed from
166
the usual FITS notation. ny is the fastest varying dimension, which is
167
usually considered the X axis in the FITS image display
168
169
*/
170
171
int stat;
172
173
if (*status > 0) return(*status);
174
175
/* H-transform */
176
stat = htrans(a, nx, ny);
177
if (stat) {
178
*status = stat;
179
return(*status);
180
}
181
182
/* digitize */
183
digitize(a, nx, ny, scale);
184
185
/* encode and write to output array */
186
187
noutmax = *nbytes; /* input value is the allocated size of the array */
188
*nbytes = 0; /* reset */
189
190
stat = encode(output, nbytes, a, nx, ny, scale);
191
192
*status = stat;
193
return(*status);
194
}
195
/* ---------------------------------------------------------------------- */
196
int _pyfits_fits_hcompress64(LONGLONG *a, int ny, int nx, int scale,
197
char *output, long *nbytes, int *status)
198
{
199
/*
200
compress the input image using the H-compress algorithm
201
202
a - input image array
203
nx - size of X axis of image
204
ny - size of Y axis of image
205
scale - quantization scale factor. Larger values results in more (lossy) compression
206
scale = 0 does lossless compression
207
output - pre-allocated array to hold the output compressed stream of bytes
208
nbyts - size of the compressed byte stream, in bytes
209
210
NOTE: the nx and ny dimensions as defined within this code are reversed from
211
the usual FITS notation. ny is the fastest varying dimension, which is
212
usually considered the X axis in the FITS image display
213
214
*/
215
216
int stat;
217
218
if (*status > 0) return(*status);
219
220
/* H-transform */
221
stat = htrans64(a, nx, ny);
222
if (stat) {
223
*status = stat;
224
return(*status);
225
}
226
227
/* digitize */
228
digitize64(a, nx, ny, scale);
229
230
/* encode and write to output array */
231
noutmax = *nbytes; /* input value is the allocated size of the array */
232
*nbytes = 0; /* reset */
233
234
stat = encode64(output, nbytes, a, nx, ny, scale);
235
236
*status = stat;
237
return(*status);
238
}
239
240
241
/* Copyright (c) 1993 Association of Universities for Research
242
* in Astronomy. All rights reserved. Produced under National
243
* Aeronautics and Space Administration Contract No. NAS5-26555.
244
*/
245
/* htrans.c H-transform of NX x NY integer image
246
*
247
* Programmer: R. White Date: 11 May 1992
248
*/
249
250
/* ######################################################################### */
251
static int htrans(int a[],int nx,int ny)
252
{
253
int nmax, log2n, h0, hx, hy, hc, nxtop, nytop, i, j, k;
254
int oddx, oddy;
255
int shift, mask, mask2, prnd, prnd2, nrnd2;
256
int s10, s00;
257
int *tmp;
258
259
/*
260
* log2n is log2 of max(nx,ny) rounded up to next power of 2
261
*/
262
nmax = (nx>ny) ? nx : ny;
263
log2n = (int) (log((float) nmax)/log(2.0)+0.5);
264
if ( nmax > (1<<log2n) ) {
265
log2n += 1;
266
}
267
/*
268
* get temporary storage for shuffling elements
269
*/
270
tmp = (int *) malloc(((nmax+1)/2)*sizeof(int));
271
if(tmp == (int *) NULL) {
272
_pyfits_ffpmsg("htrans: insufficient memory");
273
return(DATA_COMPRESSION_ERR);
274
}
275
/*
276
* set up rounding and shifting masks
277
*/
278
shift = 0;
279
mask = -2;
280
mask2 = mask << 1;
281
prnd = 1;
282
prnd2 = prnd << 1;
283
nrnd2 = prnd2 - 1;
284
/*
285
* do log2n reductions
286
*
287
* We're indexing a as a 2-D array with dimensions (nx,ny).
288
*/
289
nxtop = nx;
290
nytop = ny;
291
292
for (k = 0; k<log2n; k++) {
293
oddx = nxtop % 2;
294
oddy = nytop % 2;
295
for (i = 0; i<nxtop-oddx; i += 2) {
296
s00 = i*ny; /* s00 is index of a[i,j] */
297
s10 = s00+ny; /* s10 is index of a[i+1,j] */
298
for (j = 0; j<nytop-oddy; j += 2) {
299
/*
300
* Divide h0,hx,hy,hc by 2 (1 the first time through).
301
*/
302
h0 = (a[s10+1] + a[s10] + a[s00+1] + a[s00]) >> shift;
303
hx = (a[s10+1] + a[s10] - a[s00+1] - a[s00]) >> shift;
304
hy = (a[s10+1] - a[s10] + a[s00+1] - a[s00]) >> shift;
305
hc = (a[s10+1] - a[s10] - a[s00+1] + a[s00]) >> shift;
306
307
/*
308
* Throw away the 2 bottom bits of h0, bottom bit of hx,hy.
309
* To get rounding to be same for positive and negative
310
* numbers, nrnd2 = prnd2 - 1.
311
*/
312
a[s10+1] = hc;
313
a[s10 ] = ( (hx>=0) ? (hx+prnd) : hx ) & mask ;
314
a[s00+1] = ( (hy>=0) ? (hy+prnd) : hy ) & mask ;
315
a[s00 ] = ( (h0>=0) ? (h0+prnd2) : (h0+nrnd2) ) & mask2;
316
s00 += 2;
317
s10 += 2;
318
}
319
if (oddy) {
320
/*
321
* do last element in row if row length is odd
322
* s00+1, s10+1 are off edge
323
*/
324
h0 = (a[s10] + a[s00]) << (1-shift);
325
hx = (a[s10] - a[s00]) << (1-shift);
326
a[s10 ] = ( (hx>=0) ? (hx+prnd) : hx ) & mask ;
327
a[s00 ] = ( (h0>=0) ? (h0+prnd2) : (h0+nrnd2) ) & mask2;
328
s00 += 1;
329
s10 += 1;
330
}
331
}
332
if (oddx) {
333
/*
334
* do last row if column length is odd
335
* s10, s10+1 are off edge
336
*/
337
s00 = i*ny;
338
for (j = 0; j<nytop-oddy; j += 2) {
339
h0 = (a[s00+1] + a[s00]) << (1-shift);
340
hy = (a[s00+1] - a[s00]) << (1-shift);
341
a[s00+1] = ( (hy>=0) ? (hy+prnd) : hy ) & mask ;
342
a[s00 ] = ( (h0>=0) ? (h0+prnd2) : (h0+nrnd2) ) & mask2;
343
s00 += 2;
344
}
345
if (oddy) {
346
/*
347
* do corner element if both row and column lengths are odd
348
* s00+1, s10, s10+1 are off edge
349
*/
350
h0 = a[s00] << (2-shift);
351
a[s00 ] = ( (h0>=0) ? (h0+prnd2) : (h0+nrnd2) ) & mask2;
352
}
353
}
354
/*
355
* now shuffle in each dimension to group coefficients by order
356
*/
357
for (i = 0; i<nxtop; i++) {
358
shuffle(&a[ny*i],nytop,1,tmp);
359
}
360
for (j = 0; j<nytop; j++) {
361
shuffle(&a[j],nxtop,ny,tmp);
362
}
363
/*
364
* image size reduced by 2 (round up if odd)
365
*/
366
nxtop = (nxtop+1)>>1;
367
nytop = (nytop+1)>>1;
368
/*
369
* divisor doubles after first reduction
370
*/
371
shift = 1;
372
/*
373
* masks, rounding values double after each iteration
374
*/
375
mask = mask2;
376
prnd = prnd2;
377
mask2 = mask2 << 1;
378
prnd2 = prnd2 << 1;
379
nrnd2 = prnd2 - 1;
380
}
381
free(tmp);
382
return(0);
383
}
384
/* ######################################################################### */
385
386
static int htrans64(LONGLONG a[],int nx,int ny)
387
{
388
int nmax, log2n, nxtop, nytop, i, j, k;
389
int oddx, oddy;
390
int shift;
391
int s10, s00;
392
LONGLONG h0, hx, hy, hc, prnd, prnd2, nrnd2, mask, mask2;
393
LONGLONG *tmp;
394
395
/*
396
* log2n is log2 of max(nx,ny) rounded up to next power of 2
397
*/
398
nmax = (nx>ny) ? nx : ny;
399
log2n = (int) (log((float) nmax)/log(2.0)+0.5);
400
if ( nmax > (1<<log2n) ) {
401
log2n += 1;
402
}
403
/*
404
* get temporary storage for shuffling elements
405
*/
406
tmp = (LONGLONG *) malloc(((nmax+1)/2)*sizeof(LONGLONG));
407
if(tmp == (LONGLONG *) NULL) {
408
_pyfits_ffpmsg("htrans64: insufficient memory");
409
return(DATA_COMPRESSION_ERR);
410
}
411
/*
412
* set up rounding and shifting masks
413
*/
414
shift = 0;
415
mask = (LONGLONG) -2;
416
mask2 = mask << 1;
417
prnd = (LONGLONG) 1;
418
prnd2 = prnd << 1;
419
nrnd2 = prnd2 - 1;
420
/*
421
* do log2n reductions
422
*
423
* We're indexing a as a 2-D array with dimensions (nx,ny).
424
*/
425
nxtop = nx;
426
nytop = ny;
427
428
for (k = 0; k<log2n; k++) {
429
oddx = nxtop % 2;
430
oddy = nytop % 2;
431
for (i = 0; i<nxtop-oddx; i += 2) {
432
s00 = i*ny; /* s00 is index of a[i,j] */
433
s10 = s00+ny; /* s10 is index of a[i+1,j] */
434
for (j = 0; j<nytop-oddy; j += 2) {
435
/*
436
* Divide h0,hx,hy,hc by 2 (1 the first time through).
437
*/
438
h0 = (a[s10+1] + a[s10] + a[s00+1] + a[s00]) >> shift;
439
hx = (a[s10+1] + a[s10] - a[s00+1] - a[s00]) >> shift;
440
hy = (a[s10+1] - a[s10] + a[s00+1] - a[s00]) >> shift;
441
hc = (a[s10+1] - a[s10] - a[s00+1] + a[s00]) >> shift;
442
443
/*
444
* Throw away the 2 bottom bits of h0, bottom bit of hx,hy.
445
* To get rounding to be same for positive and negative
446
* numbers, nrnd2 = prnd2 - 1.
447
*/
448
a[s10+1] = hc;
449
a[s10 ] = ( (hx>=0) ? (hx+prnd) : hx ) & mask ;
450
a[s00+1] = ( (hy>=0) ? (hy+prnd) : hy ) & mask ;
451
a[s00 ] = ( (h0>=0) ? (h0+prnd2) : (h0+nrnd2) ) & mask2;
452
s00 += 2;
453
s10 += 2;
454
}
455
if (oddy) {
456
/*
457
* do last element in row if row length is odd
458
* s00+1, s10+1 are off edge
459
*/
460
h0 = (a[s10] + a[s00]) << (1-shift);
461
hx = (a[s10] - a[s00]) << (1-shift);
462
a[s10 ] = ( (hx>=0) ? (hx+prnd) : hx ) & mask ;
463
a[s00 ] = ( (h0>=0) ? (h0+prnd2) : (h0+nrnd2) ) & mask2;
464
s00 += 1;
465
s10 += 1;
466
}
467
}
468
if (oddx) {
469
/*
470
* do last row if column length is odd
471
* s10, s10+1 are off edge
472
*/
473
s00 = i*ny;
474
for (j = 0; j<nytop-oddy; j += 2) {
475
h0 = (a[s00+1] + a[s00]) << (1-shift);
476
hy = (a[s00+1] - a[s00]) << (1-shift);
477
a[s00+1] = ( (hy>=0) ? (hy+prnd) : hy ) & mask ;
478
a[s00 ] = ( (h0>=0) ? (h0+prnd2) : (h0+nrnd2) ) & mask2;
479
s00 += 2;
480
}
481
if (oddy) {
482
/*
483
* do corner element if both row and column lengths are odd
484
* s00+1, s10, s10+1 are off edge
485
*/
486
h0 = a[s00] << (2-shift);
487
a[s00 ] = ( (h0>=0) ? (h0+prnd2) : (h0+nrnd2) ) & mask2;
488
}
489
}
490
/*
491
* now shuffle in each dimension to group coefficients by order
492
*/
493
for (i = 0; i<nxtop; i++) {
494
shuffle64(&a[ny*i],nytop,1,tmp);
495
}
496
for (j = 0; j<nytop; j++) {
497
shuffle64(&a[j],nxtop,ny,tmp);
498
}
499
/*
500
* image size reduced by 2 (round up if odd)
501
*/
502
nxtop = (nxtop+1)>>1;
503
nytop = (nytop+1)>>1;
504
/*
505
* divisor doubles after first reduction
506
*/
507
shift = 1;
508
/*
509
* masks, rounding values double after each iteration
510
*/
511
mask = mask2;
512
prnd = prnd2;
513
mask2 = mask2 << 1;
514
prnd2 = prnd2 << 1;
515
nrnd2 = prnd2 - 1;
516
}
517
free(tmp);
518
return(0);
519
}
520
521
/* ######################################################################### */
522
static void
523
shuffle(int a[], int n, int n2, int tmp[])
524
{
525
526
/*
527
int a[]; array to shuffle
528
int n; number of elements to shuffle
529
int n2; second dimension
530
int tmp[]; scratch storage
531
*/
532
533
int i;
534
int *p1, *p2, *pt;
535
536
/*
537
* copy odd elements to tmp
538
*/
539
pt = tmp;
540
p1 = &a[n2];
541
for (i=1; i < n; i += 2) {
542
*pt = *p1;
543
pt += 1;
544
p1 += (n2+n2);
545
}
546
/*
547
* compress even elements into first half of A
548
*/
549
p1 = &a[n2];
550
p2 = &a[n2+n2];
551
for (i=2; i<n; i += 2) {
552
*p1 = *p2;
553
p1 += n2;
554
p2 += (n2+n2);
555
}
556
/*
557
* put odd elements into 2nd half
558
*/
559
pt = tmp;
560
for (i = 1; i<n; i += 2) {
561
*p1 = *pt;
562
p1 += n2;
563
pt += 1;
564
}
565
}
566
/* ######################################################################### */
567
static void
568
shuffle64(LONGLONG a[], int n, int n2, LONGLONG tmp[])
569
{
570
571
/*
572
LONGLONG a[]; array to shuffle
573
int n; number of elements to shuffle
574
int n2; second dimension
575
LONGLONG tmp[]; scratch storage
576
*/
577
578
int i;
579
LONGLONG *p1, *p2, *pt;
580
581
/*
582
* copy odd elements to tmp
583
*/
584
pt = tmp;
585
p1 = &a[n2];
586
for (i=1; i < n; i += 2) {
587
*pt = *p1;
588
pt += 1;
589
p1 += (n2+n2);
590
}
591
/*
592
* compress even elements into first half of A
593
*/
594
p1 = &a[n2];
595
p2 = &a[n2+n2];
596
for (i=2; i<n; i += 2) {
597
*p1 = *p2;
598
p1 += n2;
599
p2 += (n2+n2);
600
}
601
/*
602
* put odd elements into 2nd half
603
*/
604
pt = tmp;
605
for (i = 1; i<n; i += 2) {
606
*p1 = *pt;
607
p1 += n2;
608
pt += 1;
609
}
610
}
611
/* ######################################################################### */
612
/* ######################################################################### */
613
/* Copyright (c) 1993 Association of Universities for Research
614
* in Astronomy. All rights reserved. Produced under National
615
* Aeronautics and Space Administration Contract No. NAS5-26555.
616
*/
617
/* digitize.c digitize H-transform
618
*
619
* Programmer: R. White Date: 11 March 1991
620
*/
621
622
/* ######################################################################### */
623
static void
624
digitize(int a[], int nx, int ny, int scale)
625
{
626
int d, *p;
627
628
/*
629
* round to multiple of scale
630
*/
631
if (scale <= 1) return;
632
d=(scale+1)/2-1;
633
for (p=a; p <= &a[nx*ny-1]; p++) *p = ((*p>0) ? (*p+d) : (*p-d))/scale;
634
}
635
636
/* ######################################################################### */
637
static void
638
digitize64(LONGLONG a[], int nx, int ny, int scale)
639
{
640
LONGLONG d, *p, scale64;
641
642
/*
643
* round to multiple of scale
644
*/
645
if (scale <= 1) return;
646
d=(scale+1)/2-1;
647
scale64 = scale; /* use a 64-bit int for efficiency in the big loop */
648
649
for (p=a; p <= &a[nx*ny-1]; p++) *p = ((*p>0) ? (*p+d) : (*p-d))/scale64;
650
}
651
/* ######################################################################### */
652
/* ######################################################################### */
653
/* Copyright (c) 1993 Association of Universities for Research
654
* in Astronomy. All rights reserved. Produced under National
655
* Aeronautics and Space Administration Contract No. NAS5-26555.
656
*/
657
/* encode.c encode H-transform and write to outfile
658
*
659
* Programmer: R. White Date: 2 February 1994
660
*/
661
662
static char code_magic[2] = { (char)0xDD, (char)0x99 };
663
664
665
/* ######################################################################### */
666
static int encode(char *outfile, long *nlength, int a[], int nx, int ny, int scale)
667
{
668
669
/* FILE *outfile; - change outfile to a char array */
670
/*
671
long * nlength returned length (in bytes) of the encoded array)
672
int a[]; input H-transform array (nx,ny)
673
int nx,ny; size of H-transform array
674
int scale; scale factor for digitization
675
*/
676
int nel, nx2, ny2, i, j, k, q, vmax[3], nsign, bits_to_go;
677
unsigned char nbitplanes[3];
678
unsigned char *signbits;
679
int stat;
680
681
noutchar = 0; /* initialize the number of compressed bytes that have been written */
682
nel = nx*ny;
683
/*
684
* write magic value
685
*/
686
qwrite(outfile, code_magic, sizeof(code_magic));
687
writeint(outfile, nx); /* size of image */
688
writeint(outfile, ny);
689
writeint(outfile, scale); /* scale factor for digitization */
690
/*
691
* write first value of A (sum of all pixels -- the only value
692
* which does not compress well)
693
*/
694
writelonglong(outfile, (LONGLONG) a[0]);
695
696
a[0] = 0;
697
/*
698
* allocate array for sign bits and save values, 8 per byte
699
*/
700
signbits = (unsigned char *) malloc((nel+7)/8);
701
if (signbits == (unsigned char *) NULL) {
702
_pyfits_ffpmsg("encode: insufficient memory");
703
return(DATA_COMPRESSION_ERR);
704
}
705
nsign = 0;
706
bits_to_go = 8;
707
signbits[0] = 0;
708
for (i=0; i<nel; i++) {
709
if (a[i] > 0) {
710
/*
711
* positive element, put zero at end of buffer
712
*/
713
signbits[nsign] <<= 1;
714
bits_to_go -= 1;
715
} else if (a[i] < 0) {
716
/*
717
* negative element, shift in a one
718
*/
719
signbits[nsign] <<= 1;
720
signbits[nsign] |= 1;
721
bits_to_go -= 1;
722
/*
723
* replace a by absolute value
724
*/
725
a[i] = -a[i];
726
}
727
if (bits_to_go == 0) {
728
/*
729
* filled up this byte, go to the next one
730
*/
731
bits_to_go = 8;
732
nsign += 1;
733
signbits[nsign] = 0;
734
}
735
}
736
if (bits_to_go != 8) {
737
/*
738
* some bits in last element
739
* move bits in last byte to bottom and increment nsign
740
*/
741
signbits[nsign] <<= bits_to_go;
742
nsign += 1;
743
}
744
/*
745
* calculate number of bit planes for 3 quadrants
746
*
747
* quadrant 0=bottom left, 1=bottom right or top left, 2=top right,
748
*/
749
for (q=0; q<3; q++) {
750
vmax[q] = 0;
751
}
752
/*
753
* get maximum absolute value in each quadrant
754
*/
755
nx2 = (nx+1)/2;
756
ny2 = (ny+1)/2;
757
j=0; /* column counter */
758
k=0; /* row counter */
759
for (i=0; i<nel; i++) {
760
q = (j>=ny2) + (k>=nx2);
761
if (vmax[q] < a[i]) vmax[q] = a[i];
762
if (++j >= ny) {
763
j = 0;
764
k += 1;
765
}
766
}
767
/*
768
* now calculate number of bits for each quadrant
769
*/
770
771
/* this is a more efficient way to do this, */
772
773
774
for (q = 0; q < 3; q++) {
775
for (nbitplanes[q] = 0; vmax[q]>0; vmax[q] = vmax[q]>>1, nbitplanes[q]++) ;
776
}
777
778
779
/*
780
for (q = 0; q < 3; q++) {
781
nbitplanes[q] = (int) (log((float) (vmax[q]+1))/log(2.0)+0.5);
782
if ( (vmax[q]+1) > (1<<nbitplanes[q]) ) {
783
nbitplanes[q] += 1;
784
}
785
}
786
*/
787
788
/*
789
* write nbitplanes
790
*/
791
if (0 == qwrite(outfile, (char *) nbitplanes, sizeof(nbitplanes))) {
792
*nlength = noutchar;
793
_pyfits_ffpmsg("encode: output buffer too small");
794
return(DATA_COMPRESSION_ERR);
795
}
796
797
/*
798
* write coded array
799
*/
800
stat = doencode(outfile, a, nx, ny, nbitplanes);
801
/*
802
* write sign bits
803
*/
804
805
if (nsign > 0) {
806
807
if ( 0 == qwrite(outfile, (char *) signbits, nsign)) {
808
free(signbits);
809
*nlength = noutchar;
810
_pyfits_ffpmsg("encode: output buffer too small");
811
return(DATA_COMPRESSION_ERR);
812
}
813
}
814
815
free(signbits);
816
*nlength = noutchar;
817
818
if (noutchar >= noutmax) {
819
_pyfits_ffpmsg("encode: output buffer too small");
820
return(DATA_COMPRESSION_ERR);
821
}
822
823
return(stat);
824
}
825
/* ######################################################################### */
826
static int encode64(char *outfile, long *nlength, LONGLONG a[], int nx, int ny, int scale)
827
{
828
829
/* FILE *outfile; - change outfile to a char array */
830
/*
831
long * nlength returned length (in bytes) of the encoded array)
832
LONGLONG a[]; input H-transform array (nx,ny)
833
int nx,ny; size of H-transform array
834
int scale; scale factor for digitization
835
*/
836
int nel, nx2, ny2, i, j, k, q, nsign, bits_to_go;
837
LONGLONG vmax[3];
838
unsigned char nbitplanes[3];
839
unsigned char *signbits;
840
int stat;
841
842
noutchar = 0; /* initialize the number of compressed bytes that have been written */
843
nel = nx*ny;
844
/*
845
* write magic value
846
*/
847
qwrite(outfile, code_magic, sizeof(code_magic));
848
writeint(outfile, nx); /* size of image */
849
writeint(outfile, ny);
850
writeint(outfile, scale); /* scale factor for digitization */
851
/*
852
* write first value of A (sum of all pixels -- the only value
853
* which does not compress well)
854
*/
855
writelonglong(outfile, a[0]);
856
857
a[0] = 0;
858
/*
859
* allocate array for sign bits and save values, 8 per byte
860
*/
861
signbits = (unsigned char *) malloc((nel+7)/8);
862
if (signbits == (unsigned char *) NULL) {
863
_pyfits_ffpmsg("encode64: insufficient memory");
864
return(DATA_COMPRESSION_ERR);
865
}
866
nsign = 0;
867
bits_to_go = 8;
868
signbits[0] = 0;
869
for (i=0; i<nel; i++) {
870
if (a[i] > 0) {
871
/*
872
* positive element, put zero at end of buffer
873
*/
874
signbits[nsign] <<= 1;
875
bits_to_go -= 1;
876
} else if (a[i] < 0) {
877
/*
878
* negative element, shift in a one
879
*/
880
signbits[nsign] <<= 1;
881
signbits[nsign] |= 1;
882
bits_to_go -= 1;
883
/*
884
* replace a by absolute value
885
*/
886
a[i] = -a[i];
887
}
888
if (bits_to_go == 0) {
889
/*
890
* filled up this byte, go to the next one
891
*/
892
bits_to_go = 8;
893
nsign += 1;
894
signbits[nsign] = 0;
895
}
896
}
897
if (bits_to_go != 8) {
898
/*
899
* some bits in last element
900
* move bits in last byte to bottom and increment nsign
901
*/
902
signbits[nsign] <<= bits_to_go;
903
nsign += 1;
904
}
905
/*
906
* calculate number of bit planes for 3 quadrants
907
*
908
* quadrant 0=bottom left, 1=bottom right or top left, 2=top right,
909
*/
910
for (q=0; q<3; q++) {
911
vmax[q] = 0;
912
}
913
/*
914
* get maximum absolute value in each quadrant
915
*/
916
nx2 = (nx+1)/2;
917
ny2 = (ny+1)/2;
918
j=0; /* column counter */
919
k=0; /* row counter */
920
for (i=0; i<nel; i++) {
921
q = (j>=ny2) + (k>=nx2);
922
if (vmax[q] < a[i]) vmax[q] = a[i];
923
if (++j >= ny) {
924
j = 0;
925
k += 1;
926
}
927
}
928
/*
929
* now calculate number of bits for each quadrant
930
*/
931
932
/* this is a more efficient way to do this, */
933
934
935
for (q = 0; q < 3; q++) {
936
for (nbitplanes[q] = 0; vmax[q]>0; vmax[q] = vmax[q]>>1, nbitplanes[q]++) ;
937
}
938
939
940
/*
941
for (q = 0; q < 3; q++) {
942
nbitplanes[q] = log((float) (vmax[q]+1))/log(2.0)+0.5;
943
if ( (vmax[q]+1) > (((LONGLONG) 1)<<nbitplanes[q]) ) {
944
nbitplanes[q] += 1;
945
}
946
}
947
*/
948
949
/*
950
* write nbitplanes
951
*/
952
953
if (0 == qwrite(outfile, (char *) nbitplanes, sizeof(nbitplanes))) {
954
*nlength = noutchar;
955
_pyfits_ffpmsg("encode: output buffer too small");
956
return(DATA_COMPRESSION_ERR);
957
}
958
959
/*
960
* write coded array
961
*/
962
stat = doencode64(outfile, a, nx, ny, nbitplanes);
963
/*
964
* write sign bits
965
*/
966
967
if (nsign > 0) {
968
969
if ( 0 == qwrite(outfile, (char *) signbits, nsign)) {
970
free(signbits);
971
*nlength = noutchar;
972
_pyfits_ffpmsg("encode: output buffer too small");
973
return(DATA_COMPRESSION_ERR);
974
}
975
}
976
977
free(signbits);
978
*nlength = noutchar;
979
980
if (noutchar >= noutmax) {
981
_pyfits_ffpmsg("encode64: output buffer too small");
982
return(DATA_COMPRESSION_ERR);
983
}
984
985
return(stat);
986
}
987
/* ######################################################################### */
988
/* ######################################################################### */
989
/* Copyright (c) 1993 Association of Universities for Research
990
* in Astronomy. All rights reserved. Produced under National
991
* Aeronautics and Space Administration Contract No. NAS5-26555.
992
*/
993
/* qwrite.c Write binary data
994
*
995
* Programmer: R. White Date: 11 March 1991
996
*/
997
998
/* ######################################################################### */
999
static void
1000
writeint(char *outfile, int a)
1001
{
1002
int i;
1003
unsigned char b[4];
1004
1005
/* Write integer A one byte at a time to outfile.
1006
*
1007
* This is portable from Vax to Sun since it eliminates the
1008
* need for byte-swapping.
1009
*/
1010
for (i=3; i>=0; i--) {
1011
b[i] = a & 0x000000ff;
1012
a >>= 8;
1013
}
1014
for (i=0; i<4; i++) qwrite(outfile, (char *) &b[i],1);
1015
}
1016
1017
/* ######################################################################### */
1018
static void
1019
writelonglong(char *outfile, LONGLONG a)
1020
{
1021
int i;
1022
unsigned char b[8];
1023
1024
/* Write integer A one byte at a time to outfile.
1025
*
1026
* This is portable from Vax to Sun since it eliminates the
1027
* need for byte-swapping.
1028
*/
1029
for (i=7; i>=0; i--) {
1030
b[i] = (unsigned char) (a & 0x000000ff);
1031
a >>= 8;
1032
}
1033
for (i=0; i<8; i++) qwrite(outfile, (char *) &b[i],1);
1034
}
1035
/* ######################################################################### */
1036
static int
1037
qwrite(char *file, char buffer[], int n)
1038
{
1039
/*
1040
* write n bytes from buffer into file
1041
* returns number of bytes read (=n) if successful, <=0 if not
1042
*/
1043
1044
if (noutchar + n > noutmax) return(0); /* buffer overflow */
1045
1046
memcpy(&file[noutchar], buffer, n);
1047
noutchar += n;
1048
1049
return(n);
1050
}
1051
/* ######################################################################### */
1052
/* ######################################################################### */
1053
/* Copyright (c) 1993 Association of Universities for Research
1054
* in Astronomy. All rights reserved. Produced under National
1055
* Aeronautics and Space Administration Contract No. NAS5-26555.
1056
*/
1057
/* doencode.c Encode 2-D array and write stream of characters on outfile
1058
*
1059
* This version assumes that A is positive.
1060
*
1061
* Programmer: R. White Date: 7 May 1991
1062
*/
1063
1064
/* ######################################################################### */
1065
static int
1066
doencode(char *outfile, int a[], int nx, int ny, unsigned char nbitplanes[3])
1067
{
1068
/* char *outfile; output data stream
1069
int a[]; Array of values to encode
1070
int nx,ny; Array dimensions [nx][ny]
1071
unsigned char nbitplanes[3]; Number of bit planes in quadrants
1072
*/
1073
1074
int nx2, ny2, stat;
1075
1076
nx2 = (nx+1)/2;
1077
ny2 = (ny+1)/2;
1078
/*
1079
* Initialize bit output
1080
*/
1081
start_outputing_bits();
1082
/*
1083
* write out the bit planes for each quadrant
1084
*/
1085
stat = qtree_encode(outfile, &a[0], ny, nx2, ny2, nbitplanes[0]);
1086
1087
if (!stat)
1088
stat = qtree_encode(outfile, &a[ny2], ny, nx2, ny/2, nbitplanes[1]);
1089
1090
if (!stat)
1091
stat = qtree_encode(outfile, &a[ny*nx2], ny, nx/2, ny2, nbitplanes[1]);
1092
1093
if (!stat)
1094
stat = qtree_encode(outfile, &a[ny*nx2+ny2], ny, nx/2, ny/2, nbitplanes[2]);
1095
/*
1096
* Add zero as an EOF symbol
1097
*/
1098
output_nybble(outfile, 0);
1099
done_outputing_bits(outfile);
1100
1101
return(stat);
1102
}
1103
/* ######################################################################### */
1104
static int
1105
doencode64(char *outfile, LONGLONG a[], int nx, int ny, unsigned char nbitplanes[3])
1106
{
1107
/* char *outfile; output data stream
1108
LONGLONG a[]; Array of values to encode
1109
int nx,ny; Array dimensions [nx][ny]
1110
unsigned char nbitplanes[3]; Number of bit planes in quadrants
1111
*/
1112
1113
int nx2, ny2, stat;
1114
1115
nx2 = (nx+1)/2;
1116
ny2 = (ny+1)/2;
1117
/*
1118
* Initialize bit output
1119
*/
1120
start_outputing_bits();
1121
/*
1122
* write out the bit planes for each quadrant
1123
*/
1124
stat = qtree_encode64(outfile, &a[0], ny, nx2, ny2, nbitplanes[0]);
1125
1126
if (!stat)
1127
stat = qtree_encode64(outfile, &a[ny2], ny, nx2, ny/2, nbitplanes[1]);
1128
1129
if (!stat)
1130
stat = qtree_encode64(outfile, &a[ny*nx2], ny, nx/2, ny2, nbitplanes[1]);
1131
1132
if (!stat)
1133
stat = qtree_encode64(outfile, &a[ny*nx2+ny2], ny, nx/2, ny/2, nbitplanes[2]);
1134
/*
1135
* Add zero as an EOF symbol
1136
*/
1137
output_nybble(outfile, 0);
1138
done_outputing_bits(outfile);
1139
1140
return(stat);
1141
}
1142
/* ######################################################################### */
1143
/* ######################################################################### */
1144
/* Copyright (c) 1993 Association of Universities for Research
1145
* in Astronomy. All rights reserved. Produced under National
1146
* Aeronautics and Space Administration Contract No. NAS5-26555.
1147
*/
1148
/* BIT OUTPUT ROUTINES */
1149
1150
1151
static LONGLONG bitcount;
1152
1153
/* THE BIT BUFFER */
1154
1155
static int buffer2; /* Bits buffered for output */
1156
static int bits_to_go2; /* Number of bits free in buffer */
1157
1158
1159
/* ######################################################################### */
1160
/* INITIALIZE FOR BIT OUTPUT */
1161
1162
static void
1163
start_outputing_bits(void)
1164
{
1165
buffer2 = 0; /* Buffer is empty to start */
1166
bits_to_go2 = 8; /* with */
1167
bitcount = 0;
1168
}
1169
1170
/* ######################################################################### */
1171
/* OUTPUT N BITS (N must be <= 8) */
1172
1173
static void
1174
output_nbits(char *outfile, int bits, int n)
1175
{
1176
/* AND mask for the right-most n bits */
1177
static int mask[9] = {0, 1, 3, 7, 15, 31, 63, 127, 255};
1178
/*
1179
* insert bits at end of buffer
1180
*/
1181
buffer2 <<= n;
1182
/* buffer2 |= ( bits & ((1<<n)-1) ); */
1183
buffer2 |= ( bits & (*(mask+n)) );
1184
bits_to_go2 -= n;
1185
if (bits_to_go2 <= 0) {
1186
/*
1187
* buffer2 full, put out top 8 bits
1188
*/
1189
1190
outfile[noutchar] = ((buffer2>>(-bits_to_go2)) & 0xff);
1191
1192
if (noutchar < noutmax) noutchar++;
1193
1194
bits_to_go2 += 8;
1195
}
1196
bitcount += n;
1197
}
1198
/* ######################################################################### */
1199
/* OUTPUT a 4 bit nybble */
1200
static void
1201
output_nybble(char *outfile, int bits)
1202
{
1203
/*
1204
* insert 4 bits at end of buffer
1205
*/
1206
buffer2 = (buffer2<<4) | ( bits & 15 );
1207
bits_to_go2 -= 4;
1208
if (bits_to_go2 <= 0) {
1209
/*
1210
* buffer2 full, put out top 8 bits
1211
*/
1212
1213
outfile[noutchar] = ((buffer2>>(-bits_to_go2)) & 0xff);
1214
1215
if (noutchar < noutmax) noutchar++;
1216
1217
bits_to_go2 += 8;
1218
}
1219
bitcount += 4;
1220
}
1221
/* ############################################################################ */
1222
/* OUTPUT array of 4 BITS */
1223
1224
static void output_nnybble(char *outfile, int n, unsigned char array[])
1225
{
1226
/* pack the 4 lower bits in each element of the array into the outfile array */
1227
1228
int ii, jj, kk = 0, shift;
1229
1230
if (n == 1) {
1231
output_nybble(outfile, (int) array[0]);
1232
return;
1233
}
1234
/* forcing byte alignment doesn;t help, and even makes it go slightly slower
1235
if (bits_to_go2 != 8)
1236
output_nbits(outfile, kk, bits_to_go2);
1237
*/
1238
if (bits_to_go2 <= 4)
1239
{
1240
/* just room for 1 nybble; write it out separately */
1241
output_nybble(outfile, array[0]);
1242
kk++; /* index to next array element */
1243
1244
if (n == 2) /* only 1 more nybble to write out */
1245
{
1246
output_nybble(outfile, (int) array[1]);
1247
return;
1248
}
1249
}
1250
1251
1252
/* bits_to_go2 is now in the range 5 - 8 */
1253
shift = 8 - bits_to_go2;
1254
1255
/* now write out pairs of nybbles; this does not affect value of bits_to_go2 */
1256
jj = (n - kk) / 2;
1257
1258
if (bits_to_go2 == 8) {
1259
/* special case if nybbles are aligned on byte boundary */
1260
/* this actually seems to make very little differnece in speed */
1261
buffer2 = 0;
1262
for (ii = 0; ii < jj; ii++)
1263
{
1264
outfile[noutchar] = ((array[kk] & 15)<<4) | (array[kk+1] & 15);
1265
kk += 2;
1266
noutchar++;
1267
}
1268
} else {
1269
for (ii = 0; ii < jj; ii++)
1270
{
1271
buffer2 = (buffer2<<8) | ((array[kk] & 15)<<4) | (array[kk+1] & 15);
1272
kk += 2;
1273
1274
/*
1275
buffer2 full, put out top 8 bits
1276
*/
1277
1278
outfile[noutchar] = ((buffer2>>shift) & 0xff);
1279
noutchar++;
1280
}
1281
}
1282
1283
bitcount += (8 * (ii - 1));
1284
1285
/* write out last odd nybble, if present */
1286
if (kk != n) output_nybble(outfile, (int) array[n - 1]);
1287
1288
return;
1289
}
1290
1291
1292
/* ######################################################################### */
1293
/* FLUSH OUT THE LAST BITS */
1294
1295
static void
1296
done_outputing_bits(char *outfile)
1297
{
1298
if(bits_to_go2 < 8) {
1299
/* putc(buffer2<<bits_to_go2,outfile); */
1300
1301
outfile[noutchar] = (buffer2<<bits_to_go2);
1302
if (noutchar < noutmax) noutchar++;
1303
1304
/* count the garbage bits too */
1305
bitcount += bits_to_go2;
1306
}
1307
}
1308
/* ######################################################################### */
1309
/* ######################################################################### */
1310
/* Copyright (c) 1993 Association of Universities for Research
1311
* in Astronomy. All rights reserved. Produced under National
1312
* Aeronautics and Space Administration Contract No. NAS5-26555.
1313
*/
1314
/* qtree_encode.c Encode values in quadrant of 2-D array using binary
1315
* quadtree coding for each bit plane. Assumes array is
1316
* positive.
1317
*
1318
* Programmer: R. White Date: 15 May 1991
1319
*/
1320
1321
/*
1322
* Huffman code values and number of bits in each code
1323
*/
1324
static int code[16] =
1325
{
1326
0x3e, 0x00, 0x01, 0x08, 0x02, 0x09, 0x1a, 0x1b,
1327
0x03, 0x1c, 0x0a, 0x1d, 0x0b, 0x1e, 0x3f, 0x0c
1328
};
1329
static int ncode[16] =
1330
{
1331
6, 3, 3, 4, 3, 4, 5, 5,
1332
3, 5, 4, 5, 4, 5, 6, 4
1333
};
1334
1335
/*
1336
* variables for bit output to buffer when Huffman coding
1337
*/
1338
static int bitbuffer, bits_to_go3;
1339
1340
/*
1341
* macros to write out 4-bit nybble, Huffman code for this value
1342
*/
1343
1344
1345
/* ######################################################################### */
1346
static int
1347
qtree_encode(char *outfile, int a[], int n, int nqx, int nqy, int nbitplanes)
1348
{
1349
1350
/*
1351
int a[];
1352
int n; physical dimension of row in a
1353
int nqx; length of row
1354
int nqy; length of column (<=n)
1355
int nbitplanes; number of bit planes to output
1356
*/
1357
1358
int log2n, i, k, bit, b, bmax, nqmax, nqx2, nqy2, nx, ny;
1359
unsigned char *scratch, *buffer;
1360
1361
/*
1362
* log2n is log2 of max(nqx,nqy) rounded up to next power of 2
1363
*/
1364
nqmax = (nqx>nqy) ? nqx : nqy;
1365
log2n = (int) (log((float) nqmax)/log(2.0)+0.5);
1366
if (nqmax > (1<<log2n)) {
1367
log2n += 1;
1368
}
1369
/*
1370
* initialize buffer point, max buffer size
1371
*/
1372
nqx2 = (nqx+1)/2;
1373
nqy2 = (nqy+1)/2;
1374
bmax = (nqx2*nqy2+1)/2;
1375
/*
1376
* We're indexing A as a 2-D array with dimensions (nqx,nqy).
1377
* Scratch is 2-D with dimensions (nqx/2,nqy/2) rounded up.
1378
* Buffer is used to store string of codes for output.
1379
*/
1380
scratch = (unsigned char *) malloc(2*bmax);
1381
buffer = (unsigned char *) malloc(bmax);
1382
if ((scratch == (unsigned char *) NULL) ||
1383
(buffer == (unsigned char *) NULL)) {
1384
_pyfits_ffpmsg("qtree_encode: insufficient memory");
1385
return(DATA_COMPRESSION_ERR);
1386
}
1387
/*
1388
* now encode each bit plane, starting with the top
1389
*/
1390
for (bit=nbitplanes-1; bit >= 0; bit--) {
1391
/*
1392
* initial bit buffer
1393
*/
1394
b = 0;
1395
bitbuffer = 0;
1396
bits_to_go3 = 0;
1397
/*
1398
* on first pass copy A to scratch array
1399
*/
1400
qtree_onebit(a,n,nqx,nqy,scratch,bit);
1401
nx = (nqx+1)>>1;
1402
ny = (nqy+1)>>1;
1403
/*
1404
* copy non-zero values to output buffer, which will be written
1405
* in reverse order
1406
*/
1407
if (bufcopy(scratch,nx*ny,buffer,&b,bmax)) {
1408
/*
1409
* quadtree is expanding data,
1410
* change warning code and just fill buffer with bit-map
1411
*/
1412
write_bdirect(outfile,a,n,nqx,nqy,scratch,bit);
1413
goto bitplane_done;
1414
}
1415
/*
1416
* do log2n reductions
1417
*/
1418
for (k = 1; k<log2n; k++) {
1419
qtree_reduce(scratch,ny,nx,ny,scratch);
1420
nx = (nx+1)>>1;
1421
ny = (ny+1)>>1;
1422
if (bufcopy(scratch,nx*ny,buffer,&b,bmax)) {
1423
write_bdirect(outfile,a,n,nqx,nqy,scratch,bit);
1424
goto bitplane_done;
1425
}
1426
}
1427
/*
1428
* OK, we've got the code in buffer
1429
* Write quadtree warning code, then write buffer in reverse order
1430
*/
1431
output_nybble(outfile,0xF);
1432
if (b==0) {
1433
if (bits_to_go3>0) {
1434
/*
1435
* put out the last few bits
1436
*/
1437
output_nbits(outfile, bitbuffer & ((1<<bits_to_go3)-1),
1438
bits_to_go3);
1439
} else {
1440
/*
1441
* have to write a zero nybble if there are no 1's in array
1442
*/
1443
output_huffman(outfile,0);
1444
}
1445
} else {
1446
if (bits_to_go3>0) {
1447
/*
1448
* put out the last few bits
1449
*/
1450
output_nbits(outfile, bitbuffer & ((1<<bits_to_go3)-1),
1451
bits_to_go3);
1452
}
1453
for (i=b-1; i>=0; i--) {
1454
output_nbits(outfile,buffer[i],8);
1455
}
1456
}
1457
bitplane_done: ;
1458
}
1459
free(buffer);
1460
free(scratch);
1461
return(0);
1462
}
1463
/* ######################################################################### */
1464
static int
1465
qtree_encode64(char *outfile, LONGLONG a[], int n, int nqx, int nqy, int nbitplanes)
1466
{
1467
1468
/*
1469
LONGLONG a[];
1470
int n; physical dimension of row in a
1471
int nqx; length of row
1472
int nqy; length of column (<=n)
1473
int nbitplanes; number of bit planes to output
1474
*/
1475
1476
int log2n, i, k, bit, b, nqmax, nqx2, nqy2, nx, ny;
1477
int bmax; /* this potentially needs to be made a 64-bit int to support large arrays */
1478
unsigned char *scratch, *buffer;
1479
1480
/*
1481
* log2n is log2 of max(nqx,nqy) rounded up to next power of 2
1482
*/
1483
nqmax = (nqx>nqy) ? nqx : nqy;
1484
log2n = (int) (log((float) nqmax)/log(2.0)+0.5);
1485
if (nqmax > (1<<log2n)) {
1486
log2n += 1;
1487
}
1488
/*
1489
* initialize buffer point, max buffer size
1490
*/
1491
nqx2 = (nqx+1)/2;
1492
nqy2 = (nqy+1)/2;
1493
bmax = (( nqx2)* ( nqy2)+1)/2;
1494
/*
1495
* We're indexing A as a 2-D array with dimensions (nqx,nqy).
1496
* Scratch is 2-D with dimensions (nqx/2,nqy/2) rounded up.
1497
* Buffer is used to store string of codes for output.
1498
*/
1499
scratch = (unsigned char *) malloc(2*bmax);
1500
buffer = (unsigned char *) malloc(bmax);
1501
if ((scratch == (unsigned char *) NULL) ||
1502
(buffer == (unsigned char *) NULL)) {
1503
_pyfits_ffpmsg("qtree_encode64: insufficient memory");
1504
return(DATA_COMPRESSION_ERR);
1505
}
1506
/*
1507
* now encode each bit plane, starting with the top
1508
*/
1509
for (bit=nbitplanes-1; bit >= 0; bit--) {
1510
/*
1511
* initial bit buffer
1512
*/
1513
b = 0;
1514
bitbuffer = 0;
1515
bits_to_go3 = 0;
1516
/*
1517
* on first pass copy A to scratch array
1518
*/
1519
qtree_onebit64(a,n,nqx,nqy,scratch,bit);
1520
nx = (nqx+1)>>1;
1521
ny = (nqy+1)>>1;
1522
/*
1523
* copy non-zero values to output buffer, which will be written
1524
* in reverse order
1525
*/
1526
if (bufcopy(scratch,nx*ny,buffer,&b,bmax)) {
1527
/*
1528
* quadtree is expanding data,
1529
* change warning code and just fill buffer with bit-map
1530
*/
1531
write_bdirect64(outfile,a,n,nqx,nqy,scratch,bit);
1532
goto bitplane_done;
1533
}
1534
/*
1535
* do log2n reductions
1536
*/
1537
for (k = 1; k<log2n; k++) {
1538
qtree_reduce(scratch,ny,nx,ny,scratch);
1539
nx = (nx+1)>>1;
1540
ny = (ny+1)>>1;
1541
if (bufcopy(scratch,nx*ny,buffer,&b,bmax)) {
1542
write_bdirect64(outfile,a,n,nqx,nqy,scratch,bit);
1543
goto bitplane_done;
1544
}
1545
}
1546
/*
1547
* OK, we've got the code in buffer
1548
* Write quadtree warning code, then write buffer in reverse order
1549
*/
1550
output_nybble(outfile,0xF);
1551
if (b==0) {
1552
if (bits_to_go3>0) {
1553
/*
1554
* put out the last few bits
1555
*/
1556
output_nbits(outfile, bitbuffer & ((1<<bits_to_go3)-1),
1557
bits_to_go3);
1558
} else {
1559
/*
1560
* have to write a zero nybble if there are no 1's in array
1561
*/
1562
output_huffman(outfile,0);
1563
}
1564
} else {
1565
if (bits_to_go3>0) {
1566
/*
1567
* put out the last few bits
1568
*/
1569
output_nbits(outfile, bitbuffer & ((1<<bits_to_go3)-1),
1570
bits_to_go3);
1571
}
1572
for (i=b-1; i>=0; i--) {
1573
output_nbits(outfile,buffer[i],8);
1574
}
1575
}
1576
bitplane_done: ;
1577
}
1578
free(buffer);
1579
free(scratch);
1580
return(0);
1581
}
1582
1583
/* ######################################################################### */
1584
/*
1585
* copy non-zero codes from array to buffer
1586
*/
1587
static int
1588
bufcopy(unsigned char a[], int n, unsigned char buffer[], int *b, int bmax)
1589
{
1590
int i;
1591
1592
for (i = 0; i < n; i++) {
1593
if (a[i] != 0) {
1594
/*
1595
* add Huffman code for a[i] to buffer
1596
*/
1597
bitbuffer |= code[a[i]] << bits_to_go3;
1598
bits_to_go3 += ncode[a[i]];
1599
if (bits_to_go3 >= 8) {
1600
buffer[*b] = bitbuffer & 0xFF;
1601
*b += 1;
1602
/*
1603
* return warning code if we fill buffer
1604
*/
1605
if (*b >= bmax) return(1);
1606
bitbuffer >>= 8;
1607
bits_to_go3 -= 8;
1608
}
1609
}
1610
}
1611
return(0);
1612
}
1613
1614
/* ######################################################################### */
1615
/*
1616
* Do first quadtree reduction step on bit BIT of array A.
1617
* Results put into B.
1618
*
1619
*/
1620
static void
1621
qtree_onebit(int a[], int n, int nx, int ny, unsigned char b[], int bit)
1622
{
1623
int i, j, k;
1624
int b0, b1, b2, b3;
1625
int s10, s00;
1626
1627
/*
1628
* use selected bit to get amount to shift
1629
*/
1630
b0 = 1<<bit;
1631
b1 = b0<<1;
1632
b2 = b0<<2;
1633
b3 = b0<<3;
1634
k = 0; /* k is index of b[i/2,j/2] */
1635
for (i = 0; i<nx-1; i += 2) {
1636
s00 = n*i; /* s00 is index of a[i,j] */
1637
s10 = s00+n; /* s10 is index of a[i+1,j] */
1638
for (j = 0; j<ny-1; j += 2) {
1639
b[k] = ( ( a[s10+1] & b0)
1640
| ((a[s10 ]<<1) & b1)
1641
| ((a[s00+1]<<2) & b2)
1642
| ((a[s00 ]<<3) & b3) ) >> bit;
1643
1644
k += 1;
1645
s00 += 2;
1646
s10 += 2;
1647
}
1648
if (j < ny) {
1649
/*
1650
* row size is odd, do last element in row
1651
* s00+1,s10+1 are off edge
1652
*/
1653
b[k] = ( ((a[s10 ]<<1) & b1)
1654
| ((a[s00 ]<<3) & b3) ) >> bit;
1655
k += 1;
1656
}
1657
}
1658
if (i < nx) {
1659
/*
1660
* column size is odd, do last row
1661
* s10,s10+1 are off edge
1662
*/
1663
s00 = n*i;
1664
for (j = 0; j<ny-1; j += 2) {
1665
b[k] = ( ((a[s00+1]<<2) & b2)
1666
| ((a[s00 ]<<3) & b3) ) >> bit;
1667
k += 1;
1668
s00 += 2;
1669
}
1670
if (j < ny) {
1671
/*
1672
* both row and column size are odd, do corner element
1673
* s00+1, s10, s10+1 are off edge
1674
*/
1675
b[k] = ( ((a[s00 ]<<3) & b3) ) >> bit;
1676
k += 1;
1677
}
1678
}
1679
}
1680
/* ######################################################################### */
1681
/*
1682
* Do first quadtree reduction step on bit BIT of array A.
1683
* Results put into B.
1684
*
1685
*/
1686
static void
1687
qtree_onebit64(LONGLONG a[], int n, int nx, int ny, unsigned char b[], int bit)
1688
{
1689
int i, j, k;
1690
LONGLONG b0, b1, b2, b3;
1691
int s10, s00;
1692
1693
/*
1694
* use selected bit to get amount to shift
1695
*/
1696
b0 = ((LONGLONG) 1)<<bit;
1697
b1 = b0<<1;
1698
b2 = b0<<2;
1699
b3 = b0<<3;
1700
k = 0; /* k is index of b[i/2,j/2] */
1701
for (i = 0; i<nx-1; i += 2) {
1702
s00 = n*i; /* s00 is index of a[i,j] */
1703
s10 = s00+n; /* s10 is index of a[i+1,j] */
1704
for (j = 0; j<ny-1; j += 2) {
1705
b[k] = (unsigned char) (( ( a[s10+1] & b0)
1706
| ((a[s10 ]<<1) & b1)
1707
| ((a[s00+1]<<2) & b2)
1708
| ((a[s00 ]<<3) & b3) ) >> bit);
1709
k += 1;
1710
s00 += 2;
1711
s10 += 2;
1712
}
1713
if (j < ny) {
1714
/*
1715
* row size is odd, do last element in row
1716
* s00+1,s10+1 are off edge
1717
*/
1718
b[k] = (unsigned char) (( ((a[s10 ]<<1) & b1)
1719
| ((a[s00 ]<<3) & b3) ) >> bit);
1720
k += 1;
1721
}
1722
}
1723
if (i < nx) {
1724
/*
1725
* column size is odd, do last row
1726
* s10,s10+1 are off edge
1727
*/
1728
s00 = n*i;
1729
for (j = 0; j<ny-1; j += 2) {
1730
b[k] = (unsigned char) (( ((a[s00+1]<<2) & b2)
1731
| ((a[s00 ]<<3) & b3) ) >> bit);
1732
k += 1;
1733
s00 += 2;
1734
}
1735
if (j < ny) {
1736
/*
1737
* both row and column size are odd, do corner element
1738
* s00+1, s10, s10+1 are off edge
1739
*/
1740
b[k] = (unsigned char) (( ((a[s00 ]<<3) & b3) ) >> bit);
1741
k += 1;
1742
}
1743
}
1744
}
1745
1746
/* ######################################################################### */
1747
/*
1748
* do one quadtree reduction step on array a
1749
* results put into b (which may be the same as a)
1750
*/
1751
static void
1752
qtree_reduce(unsigned char a[], int n, int nx, int ny, unsigned char b[])
1753
{
1754
int i, j, k;
1755
int s10, s00;
1756
1757
k = 0; /* k is index of b[i/2,j/2] */
1758
for (i = 0; i<nx-1; i += 2) {
1759
s00 = n*i; /* s00 is index of a[i,j] */
1760
s10 = s00+n; /* s10 is index of a[i+1,j] */
1761
for (j = 0; j<ny-1; j += 2) {
1762
b[k] = (a[s10+1] != 0)
1763
| ( (a[s10 ] != 0) << 1)
1764
| ( (a[s00+1] != 0) << 2)
1765
| ( (a[s00 ] != 0) << 3);
1766
k += 1;
1767
s00 += 2;
1768
s10 += 2;
1769
}
1770
if (j < ny) {
1771
/*
1772
* row size is odd, do last element in row
1773
* s00+1,s10+1 are off edge
1774
*/
1775
b[k] = ( (a[s10 ] != 0) << 1)
1776
| ( (a[s00 ] != 0) << 3);
1777
k += 1;
1778
}
1779
}
1780
if (i < nx) {
1781
/*
1782
* column size is odd, do last row
1783
* s10,s10+1 are off edge
1784
*/
1785
s00 = n*i;
1786
for (j = 0; j<ny-1; j += 2) {
1787
b[k] = ( (a[s00+1] != 0) << 2)
1788
| ( (a[s00 ] != 0) << 3);
1789
k += 1;
1790
s00 += 2;
1791
}
1792
if (j < ny) {
1793
/*
1794
* both row and column size are odd, do corner element
1795
* s00+1, s10, s10+1 are off edge
1796
*/
1797
b[k] = ( (a[s00 ] != 0) << 3);
1798
k += 1;
1799
}
1800
}
1801
}
1802
1803
/* ######################################################################### */
1804
static void
1805
write_bdirect(char *outfile, int a[], int n,int nqx, int nqy, unsigned char scratch[], int bit)
1806
{
1807
1808
/*
1809
* Write the direct bitmap warning code
1810
*/
1811
output_nybble(outfile,0x0);
1812
/*
1813
* Copy A to scratch array (again!), packing 4 bits/nybble
1814
*/
1815
qtree_onebit(a,n,nqx,nqy,scratch,bit);
1816
/*
1817
* write to outfile
1818
*/
1819
/*
1820
int i;
1821
for (i = 0; i < ((nqx+1)/2) * ((nqy+1)/2); i++) {
1822
output_nybble(outfile,scratch[i]);
1823
}
1824
*/
1825
output_nnybble(outfile, ((nqx+1)/2) * ((nqy+1)/2), scratch);
1826
1827
}
1828
/* ######################################################################### */
1829
static void
1830
write_bdirect64(char *outfile, LONGLONG a[], int n,int nqx, int nqy, unsigned char scratch[], int bit)
1831
{
1832
1833
/*
1834
* Write the direct bitmap warning code
1835
*/
1836
output_nybble(outfile,0x0);
1837
/*
1838
* Copy A to scratch array (again!), packing 4 bits/nybble
1839
*/
1840
qtree_onebit64(a,n,nqx,nqy,scratch,bit);
1841
/*
1842
* write to outfile
1843
*/
1844
/*
1845
int i;
1846
for (i = 0; i < ((nqx+1)/2) * ((nqy+1)/2); i++) {
1847
output_nybble(outfile,scratch[i]);
1848
}
1849
*/
1850
output_nnybble(outfile, ((nqx+1)/2) * ((nqy+1)/2), scratch);
1851
}
Loggerhead is a web-based interface for Breezy
Version: 2.0.1