1 /* FPACK utility routines
2 SPDX-FileCopyrightText: William D. Pence <https://heasarc.gsfc.nasa.gov/fitsio/>
3 SPDX-FileCopyrightText: R. Seaman
4
5 SPDX-License-Identifier: LicenseRef-NASA-FV-License-Agreement
6 */
7
8 #include <time.h>
9 #include <float.h>
10 #include <signal.h>
11 #include <ctype.h>
12
13 /* #include "bzlib.h" only for experimental purposes */
14
15 #if defined(unix) || defined(__unix__) || defined(__unix)
16 #include <sys/time.h>
17 #endif
18
19 #include <math.h>
20 #include <fitsio.h>
21 #include <fitsio2.h>
22 #include "fpack.h"
23
24 /* these filename buffer are used to delete temporary files */
25 /* in case the program is aborted */
26 char tempfilename[SZ_STR];
27 char tempfilename2[SZ_STR];
28 char tempfilename3[SZ_STR];
29
30 /* nearest integer function */
31 # define NINT(x) ((x >= 0.) ? (int) (x + 0.5) : (int) (x - 0.5))
32 # define NSHRT(x) ((x >= 0.) ? (short) (x + 0.5) : (short) (x - 0.5))
33
34 /* define variables for measuring elapsed time */
35 clock_t scpu, ecpu;
36 long startsec; /* start of elapsed time interval */
37 int startmilli; /* start of elapsed time interval */
38
39 /* CLOCKS_PER_SEC should be defined by most compilers */
40 #if defined(CLOCKS_PER_SEC)
41 #define CLOCKTICKS CLOCKS_PER_SEC
42 #else
43 /* on SUN OS machine, CLOCKS_PER_SEC is not defined, so set its value */
44 #define CLOCKTICKS 1000000
45 #endif
46
47 FILE *outreport;
48
49 /* dimension of central image area to be sampled for test statistics */
50 int XSAMPLE = 4100;
51 int YSAMPLE = 4100;
52
53 #define UNUSED(x) (void)(x)
54 #define fp_tmpnam(suffix, rootname, tmpnam) _fp_tmpnam((char *)suffix, (char *)rootname, (char *)tmpnam)
55
56 /*--------------------------------------------------------------------------*/
fp_noop(void)57 int fp_noop (void)
58 {
59 fp_msg ("Input and output files are unchanged.\n");
60 return(0);
61 }
62 /*--------------------------------------------------------------------------*/
fp_abort_output(fitsfile * infptr,fitsfile * outfptr,int stat)63 void fp_abort_output (fitsfile *infptr, fitsfile *outfptr, int stat)
64 {
65 int status = 0, hdunum;
66 char msg[SZ_STR];
67
68 if (infptr)
69 {
70 fits_file_name(infptr, tempfilename, &status);
71 fits_get_hdu_num(infptr, &hdunum);
72
73 fits_close_file (infptr, &status);
74
75 snprintf(msg, SZ_STR,"Error processing file: %s\n", tempfilename);
76 fp_msg (msg);
77 snprintf(msg, SZ_STR," in HDU number %d\n", hdunum);
78 fp_msg (msg);
79 }
80 else
81 {
82 snprintf(msg, SZ_STR,"Error: Unable to process input file\n");
83 fp_msg(msg);
84 }
85 fits_report_error (stderr, stat);
86
87 if (outfptr) {
88 fits_delete_file(outfptr, &status);
89 fp_msg ("Input file is unchanged.\n");
90 }
91 }
92 /*--------------------------------------------------------------------------*/
fp_version(void)93 int fp_version (void)
94 {
95 float version;
96 char cfitsioversion[40];
97
98 fp_msg (FPACK_VERSION);
99 fits_get_version(&version);
100 snprintf(cfitsioversion, 40," CFITSIO version %5.3f", version);
101 fp_msg(cfitsioversion);
102 fp_msg ("\n");
103 return(0);
104 }
105 /*--------------------------------------------------------------------------*/
fp_access(char * filename)106 int fp_access (char *filename)
107 {
108 /* test if a file exists */
109
110 FILE *diskfile;
111
112 diskfile = fopen(filename, "r");
113
114 if (diskfile) {
115 fclose(diskfile);
116 return(0);
117 } else {
118 return(-1);
119 }
120 }
121 /*--------------------------------------------------------------------------*/
_fp_tmpnam(char * suffix,char * rootname,char * tmpnam)122 int _fp_tmpnam(char *suffix, char *rootname, char *tmpnam)
123 {
124 /* create temporary file name */
125
126 int maxtry = 30, ii;
127
128 if (strlen(suffix) + strlen(rootname) > SZ_STR-5) {
129 fp_msg ("Error: filename is too long to create temporary file\n"); exit (-1);
130 }
131
132 strcpy (tmpnam, rootname); /* start with rootname */
133 strcat(tmpnam, suffix); /* append the suffix */
134
135 maxtry = SZ_STR - strlen(tmpnam) - 1;
136
137 for (ii = 0; ii < maxtry; ii++) {
138 if (fp_access(tmpnam)) break; /* good, the file does not exist */
139 if (strlen(tmpnam) > SZ_STR-2)
140 {
141 fp_msg ("\nCould not create temporary file name:\n");
142 fp_msg (tmpnam);
143 fp_msg ("\n");
144 exit (-1);
145 }
146 strcat(tmpnam, "x"); /* append an x to the name, and try again */
147 }
148
149 if (ii == maxtry) {
150 fp_msg ("\nCould not create temporary file name:\n");
151 fp_msg (tmpnam);
152 fp_msg ("\n");
153 exit (-1);
154 }
155
156 return(0);
157 }
158 /*--------------------------------------------------------------------------*/
fp_init(fpstate * fpptr)159 int fp_init (fpstate *fpptr)
160 {
161 int ii;
162
163 fpptr->comptype = RICE_1;
164 fpptr->quantize_level = DEF_QLEVEL;
165 fpptr->no_dither = 0;
166 fpptr->dither_method = 1;
167 fpptr->dither_offset = 0;
168 fpptr->int_to_float = 0;
169
170 /* thresholds when using the -i2f flag */
171 fpptr->n3ratio = 2.0; /* minimum ratio of image noise sigma / q */
172 fpptr->n3min = 6.; /* minimum noise sigma. */
173
174 fpptr->scale = DEF_HCOMP_SCALE;
175 fpptr->smooth = DEF_HCOMP_SMOOTH;
176 fpptr->rescale_noise = DEF_RESCALE_NOISE;
177 fpptr->ntile[0] = (long) -1; /* -1 means extent of axis */
178
179 for (ii=1; ii < MAX_COMPRESS_DIM; ii++)
180 fpptr->ntile[ii] = (long) 1;
181
182 fpptr->to_stdout = 0;
183 fpptr->listonly = 0;
184 fpptr->clobber = 0;
185 fpptr->delete_input = 0;
186 fpptr->do_not_prompt = 0;
187 fpptr->do_checksums = 1;
188 fpptr->do_gzip_file = 0;
189 fpptr->do_tables = 0; /* this is intended for testing purposes */
190 fpptr->do_images = 1; /* can be turned off with -tableonly switch */
191 fpptr->test_all = 0;
192 fpptr->verbose = 0;
193
194 fpptr->prefix[0] = 0;
195 fpptr->extname[0] = 0;
196 fpptr->delete_suffix = 0;
197 fpptr->outfile[0] = 0;
198
199 fpptr->firstfile = 1;
200
201 /* magic number for initialization check, boolean for preflight
202 */
203 fpptr->initialized = FP_INIT_MAGIC;
204 fpptr->preflight_checked = 0;
205 return(0);
206 }
207 /*--------------------------------------------------------------------------*/
fp_list(int argc,char * argv[],fpstate fpvar)208 int fp_list (int argc, char *argv[], fpstate fpvar)
209 {
210 fitsfile *infptr;
211 char infits[SZ_STR], msg[SZ_STR];
212 int hdunum, iarg, stat=0;
213 LONGLONG sizell;
214
215 if (fpvar.initialized != FP_INIT_MAGIC) {
216 fp_msg ("Error: internal initialization error\n"); exit (-1);
217 }
218
219 for (iarg=fpvar.firstfile; iarg < argc; iarg++) {
220 strncpy (infits, argv[iarg], SZ_STR-1);
221 infits[SZ_STR-1]=0;
222
223 if (strchr (infits, '[') || strchr (infits, ']')) {
224 fp_msg ("Error: section/extension notation not supported: ");
225 fp_msg (infits); fp_msg ("\n"); exit (-1);
226 }
227
228 if (fp_access (infits) != 0) {
229 fp_msg ("Error: can't find or read input file "); fp_msg (infits);
230 fp_msg ("\n"); fp_noop (); exit (-1);
231 }
232
233 fits_open_file (&infptr, infits, READONLY, &stat);
234 if (stat) { fits_report_error (stderr, stat); exit (stat); }
235
236 /* move to the end of file, to get the total size in bytes */
237 fits_get_num_hdus (infptr, &hdunum, &stat);
238 fits_movabs_hdu (infptr, hdunum, NULL, &stat);
239 fits_get_hduaddrll(infptr, NULL, NULL, &sizell, &stat);
240
241
242 if (stat) {
243 fp_abort_output(infptr, NULL, stat);
244 }
245
246 snprintf (msg, SZ_STR,"# %s (", infits); fp_msg (msg);
247
248 #if defined(_MSC_VER)
249 /* Microsoft Visual C++ 6.0 uses '%I64d' syntax for 8-byte integers */
250 snprintf(msg, SZ_STR,"%I64d bytes)\n", sizell); fp_msg (msg);
251 #elif (USE_LL_SUFFIX == 1)
252 snprintf(msg, SZ_STR,"%lld bytes)\n", sizell); fp_msg (msg);
253 #else
254 snprintf(msg, SZ_STR,"%ld bytes)\n", sizell); fp_msg (msg);
255 #endif
256 fp_info_hdu (infptr);
257
258 fits_close_file (infptr, &stat);
259 if (stat) { fits_report_error (stderr, stat); exit (stat); }
260 }
261 return(0);
262 }
263 /*--------------------------------------------------------------------------*/
fp_info_hdu(fitsfile * infptr)264 int fp_info_hdu (fitsfile *infptr)
265 {
266 long naxes[9] = {1, 1, 1, 1, 1, 1, 1, 1, 1};
267 char msg[SZ_STR], val[SZ_CARD], com[SZ_CARD];
268 int naxis=0, hdutype, bitpix, hdupos, stat=0, ii;
269 unsigned long datasum, hdusum;
270
271 fits_movabs_hdu (infptr, 1, NULL, &stat);
272 if (stat) {
273 fp_abort_output(infptr, NULL, stat);
274 }
275
276 for (hdupos=1; ! stat; hdupos++) {
277 fits_get_hdu_type (infptr, &hdutype, &stat);
278 if (stat) {
279 fp_abort_output(infptr, NULL, stat);
280 }
281
282 /* fits_get_hdu_type calls unknown extensions "IMAGE_HDU"
283 * so consult XTENSION keyword itself
284 */
285 fits_read_keyword (infptr, "XTENSION", val, com, &stat);
286 if (stat == KEY_NO_EXIST) {
287 /* in primary HDU which by definition is an "image" */
288 stat=0; /* clear for later error handling */
289
290 } else if (stat) {
291 fp_abort_output(infptr, NULL, stat);
292
293 } else if (hdutype == IMAGE_HDU) {
294 /* that is, if XTENSION != "IMAGE" AND != "BINTABLE" */
295 if (strncmp (val+1, "IMAGE", 5) &&
296 strncmp (val+1, "BINTABLE", 5)) {
297
298 /* assign something other than any of these */
299 hdutype = IMAGE_HDU + ASCII_TBL + BINARY_TBL;
300 }
301 }
302
303 fits_get_chksum(infptr, &datasum, &hdusum, &stat);
304
305 if (hdutype == IMAGE_HDU) {
306 snprintf (msg, SZ_STR," %d IMAGE", hdupos); fp_msg (msg);
307 snprintf (msg, SZ_STR," SUMS=%lu/%lu", (unsigned long) (~((int) hdusum)), datasum); fp_msg (msg);
308
309 fits_get_img_param (infptr, 9, &bitpix, &naxis, naxes, &stat);
310
311 snprintf (msg, SZ_STR," BITPIX=%d", bitpix); fp_msg (msg);
312
313 if (naxis == 0) {
314 snprintf (msg, SZ_STR," [no_pixels]"); fp_msg (msg);
315 } else if (naxis == 1) {
316 snprintf (msg, SZ_STR," [%ld]", naxes[1]); fp_msg (msg);
317 } else {
318 snprintf (msg, SZ_STR," [%ld", naxes[0]); fp_msg (msg);
319 for (ii=1; ii < naxis; ii++) {
320 snprintf (msg, SZ_STR,"x%ld", naxes[ii]); fp_msg (msg);
321 }
322 fp_msg ("]");
323 }
324
325 if (fits_is_compressed_image (infptr, &stat)) {
326 fits_read_keyword (infptr, "ZCMPTYPE", val, com, &stat);
327
328 /* allow for quote in keyword value */
329 if (! strncmp (val+1, "RICE_1", 6))
330 fp_msg (" tiled_rice\n");
331 else if (! strncmp (val+1, "GZIP_1", 6))
332 fp_msg (" tiled_gzip_1\n");
333 else if (! strncmp (val+1, "GZIP_2", 6))
334 fp_msg (" tiled_gzip_2\n");
335 else if (! strncmp (val+1, "PLIO_1", 6))
336 fp_msg (" tiled_plio\n");
337 else if (! strncmp (val+1, "HCOMPRESS_1", 11))
338 fp_msg (" tiled_hcompress\n");
339 else
340 fp_msg (" unknown\n");
341
342 } else
343 fp_msg (" not_tiled\n");
344
345 } else if (hdutype == ASCII_TBL) {
346 snprintf (msg, SZ_STR," %d ASCII_TBL", hdupos); fp_msg (msg);
347 snprintf (msg, SZ_STR," SUMS=%lu/%lu\n", (unsigned long) (~((int) hdusum)), datasum); fp_msg (msg);
348
349 } else if (hdutype == BINARY_TBL) {
350 snprintf (msg, SZ_STR," %d BINARY_TBL", hdupos); fp_msg (msg);
351 snprintf (msg, SZ_STR," SUMS=%lu/%lu\n", (unsigned long) (~((int) hdusum)), datasum); fp_msg (msg);
352
353 } else {
354 snprintf (msg, SZ_STR," %d OTHER", hdupos); fp_msg (msg);
355 snprintf (msg, SZ_STR," SUMS=%lu/%lu", (unsigned long) (~((int) hdusum)), datasum); fp_msg (msg);
356 snprintf (msg, SZ_STR," %s\n", val); fp_msg (msg);
357 }
358
359 fits_movrel_hdu (infptr, 1, NULL, &stat);
360 }
361 return(0);
362 }
363
364 /*--------------------------------------------------------------------------*/
fp_preflight(int argc,char * argv[],int unpack,fpstate * fpptr)365 int fp_preflight (int argc, char *argv[], int unpack, fpstate *fpptr)
366 {
367 char infits[SZ_STR], outfits[SZ_STR];
368 int iarg, namelen, nfiles = 0;
369
370 if (fpptr->initialized != FP_INIT_MAGIC) {
371 fp_msg ("Error: internal initialization error\n"); exit (-1);
372 }
373
374 for (iarg=fpptr->firstfile; iarg < argc; iarg++) {
375
376 outfits[0] = '\0';
377
378 if (strlen(argv[iarg]) > SZ_STR - 4) { /* allow for .fz or .gz suffix */
379 fp_msg ("Error: input file name\n "); fp_msg (argv[iarg]);
380 fp_msg ("\n is too long\n"); fp_noop (); exit (-1);
381 }
382
383 strncpy (infits, argv[iarg], SZ_STR-1);
384 if (infits[0] == '-' && infits[1] != '\0') {
385 /* don't interpret this as intending to read input file from stdin */
386 fp_msg ("Error: invalid input file name\n "); fp_msg (argv[iarg]);
387 fp_msg ("\n"); fp_noop (); exit (-1);
388 }
389
390 if (strchr (infits, '[') || strchr (infits, ']')) {
391 fp_msg ("Error: section/extension notation not supported: ");
392 fp_msg (infits); fp_msg ("\n"); fp_noop (); exit (-1);
393 }
394
395 if (unpack) {
396 /* ********** This section applies to funpack ************ */
397
398 /* check that input file exists */
399 if (infits[0] != '-') { /* if not reading from stdin stream */
400 if (fp_access (infits) != 0) { /* if not, then check if */
401 strcat(infits, ".fz"); /* a .fz version exsits */
402 if (fp_access (infits) != 0) {
403 namelen = strlen(infits);
404 infits[namelen - 3] = '\0'; /* remove the .fz suffix */
405 fp_msg ("Error: can't find or read input file "); fp_msg (infits);
406 fp_msg ("\n"); fp_noop (); exit (-1);
407 }
408 } else { /* make sure a .fz version of the same file doesn't exist */
409 namelen = strlen(infits);
410 strcat(infits, ".fz");
411 if (fp_access (infits) == 0) {
412 infits[namelen] = '\0'; /* remove the .fz suffix */
413 fp_msg ("Error: ambiguous input file name. Which file should be unpacked?:\n ");
414 fp_msg (infits); fp_msg ("\n ");
415 fp_msg (infits); fp_msg (".fz\n");
416 fp_noop (); exit (-1);
417 } else {
418 infits[namelen] = '\0'; /* remove the .fz suffix */
419 }
420 }
421 }
422
423 /* if writing to stdout, then we are all done */
424 if (fpptr->to_stdout) {
425 continue;
426 }
427
428 if (fpptr->outfile[0]) { /* user specified output file name */
429 nfiles++;
430 if (nfiles > 1) {
431 fp_msg ("Error: cannot use same output file name for multiple files:\n ");
432 fp_msg (fpptr->outfile);
433 fp_msg ("\n"); fp_noop (); exit (-1);
434 }
435
436 /* check that output file doesn't exist */
437 if (fp_access (fpptr->outfile) == 0) {
438 fp_msg ("Error: output file already exists:\n ");
439 fp_msg (fpptr->outfile);
440 fp_msg ("\n "); fp_noop (); exit (-1);
441 }
442 continue;
443 }
444
445 /* construct output file name to test */
446 if (fpptr->prefix[0]) {
447 if (strlen(fpptr->prefix) + strlen(infits) > SZ_STR - 1) {
448 fp_msg ("Error: output file name for\n "); fp_msg (infits);
449 fp_msg ("\n is too long with the prefix\n"); fp_noop (); exit (-1);
450 }
451 strcat(outfits,fpptr->prefix);
452 }
453
454 /* construct output file name */
455 if (infits[0] == '-') {
456 strcpy(outfits, "output.fits");
457 } else {
458 strcpy(outfits, infits);
459 }
460
461 /* remove .gz suffix, if present (output is not gzipped) */
462 namelen = strlen(outfits);
463 if ( !strcmp(".gz", outfits + namelen - 3) ) {
464 outfits[namelen - 3] = '\0';
465 }
466
467 /* check for .fz suffix that is sometimes required */
468 /* and remove it if present */
469 if (infits[0] != '-') { /* if not reading from stdin stream */
470 namelen = strlen(outfits);
471 if ( !strcmp(".fz", outfits + namelen - 3) ) { /* suffix is present */
472 outfits[namelen - 3] = '\0';
473 } else if (fpptr->delete_suffix) { /* required suffix is missing */
474 fp_msg ("Error: input compressed file "); fp_msg (infits);
475 fp_msg ("\n does not have the default .fz suffix.\n");
476 fp_noop (); exit (-1);
477 }
478 }
479
480 /* if infits != outfits, make sure outfits doesn't already exist */
481 if (strcmp(infits, outfits)) {
482 if (fp_access (outfits) == 0) {
483 fp_msg ("Error: output file already exists:\n "); fp_msg (outfits);
484 fp_msg ("\n "); fp_noop (); exit (-1);
485 }
486 }
487
488 /* if gzipping the output, make sure .gz file doesn't exist */
489 if (fpptr->do_gzip_file) {
490 if (strlen(outfits)+3 > SZ_STR-1)
491 {
492 fp_msg ("Error: output file name too long:\n "); fp_msg (outfits);
493 fp_msg ("\n "); fp_noop (); exit (-1);
494 }
495 strcat(outfits, ".gz");
496 if (fp_access (outfits) == 0) {
497 fp_msg ("Error: output file already exists:\n "); fp_msg (outfits);
498 fp_msg ("\n "); fp_noop (); exit (-1);
499 }
500 namelen = strlen(outfits);
501 outfits[namelen - 3] = '\0'; /* remove the .gz suffix again */
502 }
503 } else {
504 /* ********** This section applies to fpack ************ */
505
506 /* check that input file exists */
507 if (infits[0] != '-') { /* if not reading from stdin stream */
508 if (fp_access (infits) != 0) { /* if not, then check if */
509 if (strlen(infits)+3 > SZ_STR-1)
510 {
511 fp_msg ("Error: input file name too long:\n "); fp_msg (infits);
512 fp_msg ("\n "); fp_noop (); exit (-1);
513 }
514 strcat(infits, ".gz"); /* a gzipped version exsits */
515 if (fp_access (infits) != 0) {
516 namelen = strlen(infits);
517 infits[namelen - 3] = '\0'; /* remove the .gz suffix */
518 fp_msg ("Error: can't find or read input file "); fp_msg (infits);
519 fp_msg ("\n"); fp_noop (); exit (-1);
520 }
521 }
522 }
523
524 /* make sure the file to pack does not already have a .fz suffix */
525 namelen = strlen(infits);
526 if ( !strcmp(".fz", infits + namelen - 3) ) {
527 fp_msg ("Error: fpack input file already has '.fz' suffix\n" ); fp_msg (infits);
528 fp_msg ("\n"); fp_noop (); exit (-1);
529 }
530
531 /* if writing to stdout, or just testing the files, then we are all done */
532 if (fpptr->to_stdout || fpptr->test_all) {
533 continue;
534 }
535
536 /* construct output file name */
537 if (infits[0] == '-') {
538 strcpy(outfits, "input.fits");
539 } else {
540 strcpy(outfits, infits);
541 }
542
543 /* remove .gz suffix, if present (output is not gzipped) */
544 namelen = strlen(outfits);
545 if ( !strcmp(".gz", outfits + namelen - 3) ) {
546 outfits[namelen - 3] = '\0';
547 }
548
549 /* remove .imh suffix (IRAF format image), and replace with .fits */
550 namelen = strlen(outfits);
551 if ( !strcmp(".imh", outfits + namelen - 4) ) {
552 outfits[namelen - 4] = '\0';
553 strcat(outfits, ".fits");
554 }
555
556 /* If not clobbering the input file, add .fz suffix to output name */
557 if (! fpptr->clobber)
558 strcat(outfits, ".fz");
559
560 /* if infits != outfits, make sure outfits doesn't already exist */
561 if (strcmp(infits, outfits)) {
562 if (fp_access (outfits) == 0) {
563 fp_msg ("Error: output file already exists:\n "); fp_msg (outfits);
564 fp_msg ("\n "); fp_noop (); exit (-1);
565 }
566 }
567 } /* end of fpack section */
568 }
569
570 fpptr->preflight_checked++;
571 return(0);
572 }
573
574 /*--------------------------------------------------------------------------*/
575 /* must run fp_preflight() before fp_loop()
576 */
fp_loop(int argc,char * argv[],int unpack,char * output_filename,fpstate fpvar)577 int fp_loop (int argc, char *argv[], int unpack, char *output_filename, fpstate fpvar)
578 {
579 char infits[SZ_STR], outfits[SZ_STR];
580 char temp[SZ_STR], answer[30];
581 char valchar[]="0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz.#()+,-_@[]/^{}";
582 int ichar=0, outlen=0, iarg, islossless, namelen, iraf_infile = 0, status = 0, ifail;
583
584 if (fpvar.initialized != FP_INIT_MAGIC) {
585 fp_msg ("Error: internal initialization error\n"); exit (-1);
586 } else if (! fpvar.preflight_checked) {
587 fp_msg ("Error: internal preflight error\n"); exit (-1);
588 }
589
590 if (fpvar.test_all && fpvar.outfile[0]) {
591 outreport = fopen(fpvar.outfile, "w");
592 fprintf(outreport," Filename Extension BITPIX NAXIS1 NAXIS2 Size N_nulls Minval Maxval Mean Sigm Noise1 Noise2 Noise3 Noise5 T_whole T_rowbyrow ");
593 fprintf(outreport,"[Comp_ratio, Pack_cpu, Unpack_cpu, Lossless readtimes] (repeated for Rice, Hcompress, and GZIP)\n");
594 }
595
596
597 tempfilename[0] = '\0';
598 tempfilename2[0] = '\0';
599 tempfilename3[0] = '\0';
600
601 /* set up signal handler to delete temporary file on abort */
602 #ifdef SIGINT
603 if (signal(SIGINT, SIG_IGN) != SIG_IGN) {
604 (void) signal(SIGINT, abort_fpack);
605 }
606 #endif
607
608 #ifdef SIGTERM
609 if (signal(SIGTERM, SIG_IGN) != SIG_IGN) {
610 (void) signal(SIGTERM, abort_fpack);
611 }
612 #endif
613
614 #ifdef SIGHUP
615 if (signal(SIGHUP, SIG_IGN) != SIG_IGN) {
616 (void) signal(SIGHUP, abort_fpack);
617 }
618 #endif
619
620 for (iarg=fpvar.firstfile; iarg < argc; iarg++) {
621
622 temp[0] = '\0';
623 outfits[0] = '\0';
624 islossless = 1;
625
626 strncpy (infits, argv[iarg], SZ_STR - 1);
627 infits[SZ_STR-1]=0;
628
629 if (unpack) {
630 /* ********** This section applies to funpack ************ */
631
632 /* find input file */
633 if (infits[0] != '-') { /* if not reading from stdin stream */
634 if (fp_access (infits) != 0) { /* if not, then */
635 strcat(infits, ".fz"); /* a .fz version must exsit */
636 }
637 }
638
639 if (fpvar.to_stdout) {
640 strcpy(outfits, "-");
641
642 } else if (fpvar.outfile[0]) { /* user specified output file name */
643 strcpy(outfits, fpvar.outfile);
644
645 } else {
646 /* construct output file name */
647 if (fpvar.prefix[0]) {
648 strcat(outfits,fpvar.prefix);
649 }
650
651 /* construct output file name */
652 if (infits[0] == '-') {
653 strcpy(outfits, "output.fits");
654 } else {
655 /*strcpy(outfits, infits);*/
656 strcpy(outfits, output_filename);
657 }
658
659 /* remove .gz suffix, if present (output is not gzipped) */
660 namelen = strlen(outfits);
661 if ( !strcmp(".gz", outfits + namelen - 3) ) {
662 outfits[namelen - 3] = '\0';
663 }
664
665 /* check for .fz suffix that is sometimes required */
666 /* and remove it if present */
667 namelen = strlen(outfits);
668 if ( !strcmp(".fz", outfits + namelen - 3) ) { /* suffix is present */
669 outfits[namelen - 3] = '\0';
670 }
671 }
672
673 } else {
674 /* ********** This section applies to fpack ************ */
675
676 if (fpvar.to_stdout) {
677 strcpy(outfits, "-");
678 } else if (! fpvar.test_all) {
679
680 /* construct output file name */
681 if (infits[0] == '-') {
682 strcpy(outfits, "input.fits");
683 } else {
684 strcpy(outfits, output_filename);
685 }
686
687 /* remove .gz suffix, if present (output is not gzipped) */
688 namelen = strlen(outfits);
689 if ( !strcmp(".gz", outfits + namelen - 3) ) {
690 outfits[namelen - 3] = '\0';
691 }
692
693 /* remove .imh suffix (IRAF format image), and replace with .fits */
694 namelen = strlen(outfits);
695 if ( !strcmp(".imh", outfits + namelen - 4) ) {
696 outfits[namelen - 4] = '\0';
697 strcat(outfits, ".fits");
698 iraf_infile = 1; /* this is an IRAF format input file */
699 /* change the output name to "NAME.fits.fz" */
700 }
701
702 /* If not clobbering the input file, add .fz suffix to output name */
703 if (! fpvar.clobber)
704 strcat(outfits, ".fz");
705 }
706 }
707
708 strncpy(temp, outfits, SZ_STR-1);
709 temp[SZ_STR-1]=0;
710
711 if (infits[0] != '-') { /* if not reading from stdin stream */
712 if (!strcmp(infits, outfits) ) { /* are input and output names the same? */
713
714 /* clobber the input file with the output file with the same name */
715 if (! fpvar.clobber) {
716 fp_msg ("\nError: must use -F flag to clobber input file.\n");
717 exit (-1);
718 }
719
720 /* create temporary file name in the output directory (same as input directory)*/
721 fp_tmpnam("Tmp1", infits, outfits);
722
723 strcpy(tempfilename, outfits); /* store temp file name, in case of abort */
724 }
725 }
726
727
728 /* *************** now do the real work ********************* */
729
730 if (fpvar.verbose && ! fpvar.to_stdout)
731 printf("%s ", infits);
732
733 if (fpvar.test_all) { /* compare all the algorithms */
734
735 /* create 2 temporary file names, in the CWD */
736 fp_tmpnam("Tmpfile1", "", tempfilename);
737 fp_tmpnam("Tmpfile2", "", tempfilename2);
738
739 fp_test (infits, tempfilename, tempfilename2, fpvar);
740
741 remove(tempfilename);
742 tempfilename[0] = '\0'; /* clear the temp file name */
743 remove(tempfilename2);
744 tempfilename2[0] = '\0';
745 continue;
746
747 } else if (unpack) {
748 if (fpvar.to_stdout) {
749 /* unpack the input file to the stdout stream */
750 fp_unpack (infits, outfits, fpvar);
751 } else {
752 /* unpack to temporary file, so other tasks can't open it until it is renamed */
753
754 /* create temporary file name, in the output directory */
755 fp_tmpnam("Tmp2", outfits, tempfilename2);
756
757 /* unpack the input file to the temporary file */
758 fp_unpack (infits, tempfilename2, fpvar);
759
760 /* rename the temporary file to it's real name */
761 ifail = rename(tempfilename2, outfits);
762 if (ifail) {
763 fp_msg("Failed to rename temporary file name:\n ");
764 fp_msg(tempfilename2);
765 fp_msg(" -> ");
766 fp_msg(outfits);
767 fp_msg("\n");
768 exit (-1);
769 } else {
770 tempfilename2[0] = '\0'; /* clear temporary file name */
771 }
772 }
773 } else {
774 fp_pack (infits, outfits, fpvar, &islossless);
775 }
776
777 if (fpvar.to_stdout) {
778 continue;
779 }
780
781 /* ********** clobber and/or delete files, if needed ************** */
782
783 if (!strcmp(infits, temp) && fpvar.clobber ) {
784
785 if (!islossless && ! fpvar.do_not_prompt) {
786 fp_msg ("\nFile ");
787 fp_msg (infits);
788 fp_msg ("\nwas compressed with a LOSSY method. Overwrite the\n");
789 fp_msg ("original file with the compressed version? (Y/N) ");
790 if (fgets(answer, 29, stdin) && answer[0] != 'Y' && answer[0] != 'y') {
791 fp_msg ("\noriginal file NOT overwritten!\n");
792 remove(outfits);
793 continue;
794 }
795 }
796
797 if (iraf_infile) { /* special case of deleting an IRAF format header and pixel file */
798 if (fits_delete_iraf_file(infits, &status)) {
799 fp_msg("\nError deleting IRAF .imh and .pix files.\n");
800 fp_msg(infits); fp_msg ("\n"); exit (-1);
801 }
802 }
803
804 #if defined(unix) || defined(__unix__) || defined(__unix)
805 /* rename clobbers input on Unix platforms */
806 if (rename (outfits, temp) != 0) {
807 fp_msg ("\nError renaming tmp file to ");
808 fp_msg (temp); fp_msg ("\n"); exit (-1);
809 }
810 #else
811 /* rename DOES NOT clobber existing files on Windows platforms */
812 /* so explicitly remove any existing file before renaming the file */
813 remove(temp);
814 if (rename (outfits, temp) != 0) {
815 fp_msg ("\nError renaming tmp file to ");
816 fp_msg (temp); fp_msg ("\n"); exit (-1);
817 }
818 #endif
819
820 tempfilename[0] = '\0'; /* clear temporary file name */
821 strcpy(outfits, temp);
822
823 } else if (fpvar.clobber || fpvar.delete_input) { /* delete the input file */
824 if (!islossless && !fpvar.do_not_prompt) { /* user did not turn off delete prompt */
825 fp_msg ("\nFile ");
826 fp_msg (infits);
827 fp_msg ("\nwas compressed with a LOSSY method. \n");
828 fp_msg ("Delete the original file? (Y/N) ");
829 if (fgets(answer, 29, stdin) && answer[0] != 'Y' && answer[0] != 'y') { /* user abort */
830 fp_msg ("\noriginal file NOT deleted!\n");
831 } else {
832 if (iraf_infile) { /* special case of deleting an IRAF format header and pixel file */
833 if (fits_delete_iraf_file(infits, &status)) {
834 fp_msg("\nError deleting IRAF .imh and .pix files.\n");
835 fp_msg(infits); fp_msg ("\n"); exit (-1);
836 }
837 } else if (remove(infits) != 0) { /* normal case of deleting input FITS file */
838 fp_msg ("\nError deleting input file ");
839 fp_msg (infits); fp_msg ("\n"); exit (-1);
840 }
841 }
842 } else { /* user said don't prompt, so just delete the input file */
843 if (iraf_infile) { /* special case of deleting an IRAF format header and pixel file */
844 if (fits_delete_iraf_file(infits, &status)) {
845 fp_msg("\nError deleting IRAF .imh and .pix files.\n");
846 fp_msg(infits); fp_msg ("\n"); exit (-1);
847 }
848 } else if (remove(infits) != 0) { /* normal case of deleting input FITS file */
849 fp_msg ("\nError deleting input file ");
850 fp_msg (infits); fp_msg ("\n"); exit (-1);
851 }
852 }
853 }
854 iraf_infile = 0;
855
856 if (fpvar.do_gzip_file) { /* gzip the output file */
857 strcpy(temp, "gzip -1 ");
858 outlen = strlen(outfits);
859 if (outlen + 8 > SZ_STR-1)
860 {
861 fp_msg("\nError: Output file name is too long.\n");
862 exit(-1);
863 }
864 for (ichar=0; ichar < outlen; ++ichar)
865 {
866 if (!strchr(valchar, outfits[ichar]))
867 {
868 fp_msg("\n Error: Invalid characters in output file name.\n");
869 exit(-1);
870 }
871 }
872 strcat(temp,outfits);
873 int rc = system(temp);
874 UNUSED(rc);
875 strcat(outfits, ".gz"); /* only possibible with funpack */
876 }
877
878 if (fpvar.verbose && ! fpvar.to_stdout)
879 printf("-> %s\n", outfits);
880
881 }
882
883 if (fpvar.test_all && fpvar.outfile[0])
884 fclose(outreport);
885 return(0);
886 }
887
888 /*--------------------------------------------------------------------------*/
889 /* fp_pack assumes the output file does not exist (checked by preflight)
890 */
fp_pack(char * infits,char * outfits,fpstate fpvar,int * islossless)891 int fp_pack (char *infits, char *outfits, fpstate fpvar, int *islossless)
892 {
893 fitsfile *infptr, *outfptr;
894 int stat=0;
895
896 fits_open_file (&infptr, infits, READONLY, &stat);
897 if (stat)
898 {
899 fits_report_error (stderr, stat);
900 return -1;
901 }
902
903 fits_create_file (&outfptr, outfits, &stat);
904 if (stat)
905 {
906 fp_abort_output(infptr, outfptr, stat);
907 return -1;
908 }
909
910 while (! stat) {
911
912 /* LOOP OVER EACH HDU */
913
914 fits_set_lossy_int (outfptr, fpvar.int_to_float, &stat);
915 fits_set_compression_type (outfptr, fpvar.comptype, &stat);
916 fits_set_tile_dim (outfptr, 6, fpvar.ntile, &stat);
917
918 if (fpvar.no_dither)
919 fits_set_quantize_method(outfptr, -1, &stat);
920 else
921 fits_set_quantize_method(outfptr, fpvar.dither_method, &stat);
922
923 fits_set_quantize_level (outfptr, fpvar.quantize_level, &stat);
924 fits_set_dither_offset(outfptr, fpvar.dither_offset, &stat);
925 fits_set_hcomp_scale (outfptr, fpvar.scale, &stat);
926 fits_set_hcomp_smooth (outfptr, fpvar.smooth, &stat);
927
928 fp_pack_hdu (infptr, outfptr, fpvar, islossless, &stat);
929
930 if (fpvar.do_checksums) {
931 fits_write_chksum (outfptr, &stat);
932 }
933
934 fits_movrel_hdu (infptr, 1, NULL, &stat);
935 }
936
937 if (stat == END_OF_FILE) stat = 0;
938
939 /* set checksum for case of newly created primary HDU */
940
941 if (fpvar.do_checksums) {
942 fits_movabs_hdu (outfptr, 1, NULL, &stat);
943 fits_write_chksum (outfptr, &stat);
944 }
945
946 if (stat) {
947 fp_abort_output(infptr, outfptr, stat);
948 }
949
950 fits_close_file (outfptr, &stat);
951 fits_close_file (infptr, &stat);
952
953 return(0);
954 }
955
956 /*--------------------------------------------------------------------------*/
957 /* fp_unpack assumes the output file does not exist
958 */
fp_unpack(char * infits,char * outfits,fpstate fpvar)959 int fp_unpack (char *infits, char *outfits, fpstate fpvar)
960 {
961 fitsfile *infptr, *outfptr;
962 int stat=0, hdutype, extnum, single = 0;
963 char *loc, *hduloc = 0, hduname[SZ_STR] = { 0 };
964
965 fits_open_file (&infptr, infits, READONLY, &stat);
966 fits_create_file (&outfptr, outfits, &stat);
967
968 if (stat) {
969 fp_abort_output(infptr, outfptr, stat);
970 }
971
972 if (fpvar.extname[0]) { /* unpack a list of HDUs? */
973
974 /* move to the first HDU in the list */
975 hduloc = fpvar.extname;
976 loc = strchr(hduloc, ','); /* look for 'comma' delimiter between names */
977
978 if (loc)
979 *loc = '\0'; /* terminate the first name in the string */
980
981 strcpy(hduname, hduloc); /* copy the first name into temporary string */
982
983 if (loc)
984 hduloc = loc + 1; /* advance to the beginning of the next name, if any */
985 else {
986 hduloc += strlen(hduname); /* end of the list */
987 single = 1; /* only 1 HDU is being unpacked */
988 }
989
990 if (isdigit( (int) hduname[0]) ) {
991 extnum = strtol(hduname, &loc, 10); /* read the string as an integer */
992
993 /* check for junk following the integer */
994 if (*loc == '\0' ) /* no junk, so move to this HDU number (+1) */
995 {
996 fits_movabs_hdu(infptr, extnum + 1, &hdutype, &stat); /* move to HDU number */
997 if (hdutype != IMAGE_HDU)
998 stat = NOT_IMAGE;
999
1000 } else { /* the string is not an integer, so must be the column name */
1001 hdutype = IMAGE_HDU;
1002 fits_movnam_hdu(infptr, hdutype, hduname, 0, &stat);
1003 }
1004 }
1005 else
1006 {
1007 /* move to the named image extension */
1008 hdutype = IMAGE_HDU;
1009 fits_movnam_hdu(infptr, hdutype, hduname, 0, &stat);
1010 }
1011 }
1012
1013 if (stat) {
1014 fp_msg ("Unable to find and move to extension '");
1015 fp_msg(hduname);
1016 fp_msg("'\n");
1017 fp_abort_output(infptr, outfptr, stat);
1018 }
1019
1020 while (! stat) {
1021
1022 if (single)
1023 stat = -1; /* special status flag to force output primary array */
1024
1025 fp_unpack_hdu (infptr, outfptr, fpvar, &stat);
1026
1027 if (fpvar.do_checksums) {
1028 fits_write_chksum (outfptr, &stat);
1029 }
1030
1031 /* move to the next HDU */
1032 if (fpvar.extname[0]) { /* unpack a list of HDUs? */
1033
1034 if (!(*hduloc)) {
1035 stat = END_OF_FILE; /* we reached the end of the list */
1036 } else {
1037 /* parse the next HDU name and move to it */
1038 loc = strchr(hduloc, ',');
1039
1040 if (loc) /* look for 'comma' delimiter between names */
1041 *loc = '\0'; /* terminate the first name in the string */
1042
1043 strcpy(hduname, hduloc); /* copy the next name into temporary string */
1044
1045 if (loc)
1046 hduloc = loc + 1; /* advance to the beginning of the next name, if any */
1047 else
1048 *hduloc = '\0'; /* end of the list */
1049
1050 if (isdigit( (int) hduname[0]) ) {
1051 extnum = strtol(hduname, &loc, 10); /* read the string as an integer */
1052
1053 /* check for junk following the integer */
1054 if (*loc == '\0' ) /* no junk, so move to this HDU number (+1) */
1055 {
1056 fits_movabs_hdu(infptr, extnum + 1, &hdutype, &stat); /* move to HDU number */
1057 if (hdutype != IMAGE_HDU)
1058 stat = NOT_IMAGE;
1059
1060 } else { /* the string is not an integer, so must be the column name */
1061 hdutype = IMAGE_HDU;
1062 fits_movnam_hdu(infptr, hdutype, hduname, 0, &stat);
1063 }
1064
1065 } else {
1066 /* move to the named image extension */
1067 hdutype = IMAGE_HDU;
1068 fits_movnam_hdu(infptr, hdutype, hduname, 0, &stat);
1069 }
1070
1071 if (stat) {
1072 fp_msg ("Unable to find and move to extension '");
1073 fp_msg(hduname);
1074 fp_msg("'\n");
1075 }
1076 }
1077 } else {
1078 /* increment to the next HDU */
1079 fits_movrel_hdu (infptr, 1, NULL, &stat);
1080 }
1081 }
1082
1083 if (stat == END_OF_FILE) stat = 0;
1084
1085 /* set checksum for case of newly created primary HDU
1086 */
1087 if (fpvar.do_checksums) {
1088 fits_movabs_hdu (outfptr, 1, NULL, &stat);
1089 fits_write_chksum (outfptr, &stat);
1090 }
1091
1092
1093 if (stat) {
1094 fp_abort_output(infptr, outfptr, stat);
1095 }
1096
1097 fits_close_file (outfptr, &stat);
1098 fits_close_file (infptr, &stat);
1099
1100 return(0);
1101 }
1102
1103 /*--------------------------------------------------------------------------*/
1104 /* fp_test assumes the output files do not exist
1105 */
fp_test(char * infits,char * outfits,char * outfits2,fpstate fpvar)1106 int fp_test (char *infits, char *outfits, char *outfits2, fpstate fpvar)
1107 {
1108 fitsfile *inputfptr, *infptr, *outfptr, *outfptr2, *tempfile;
1109
1110 long naxes[9] = {1, 1, 1, 1, 1, 1, 1, 1, 1};
1111 int stat=0, totpix=0, naxis=0, ii, hdutype, bitpix = 0, extnum = 0, len;
1112 int tstatus = 0, hdunum, rescale_flag, bpix = 8, ncols;
1113 char dtype[8], dimen[100];
1114 double bscale, rescale, noisemin;
1115 long headstart, datastart, dataend;
1116 float origdata = 0., whole_cpu, whole_elapse, row_elapse, row_cpu, xbits;
1117
1118 LONGLONG nrows;
1119 /* structure to hold image statistics (defined in fpack.h) */
1120 imgstats imagestats = { 0 };
1121
1122 fits_open_file (&inputfptr, infits, READONLY, &stat);
1123 fits_create_file (&outfptr, outfits, &stat);
1124 fits_create_file (&outfptr2, outfits2, &stat);
1125
1126 if (stat) { fits_report_error (stderr, stat); exit (stat); }
1127
1128 while (! stat) {
1129
1130 /* LOOP OVER EACH HDU */
1131 rescale_flag = 0;
1132 fits_get_hdu_type (inputfptr, &hdutype, &stat);
1133
1134 if (hdutype == IMAGE_HDU) {
1135 fits_get_img_param (inputfptr, 9, &bitpix, &naxis, naxes, &stat);
1136 for (totpix=1, ii=0; ii < 9; ii++) totpix *= naxes[ii];
1137 }
1138
1139 if (!fits_is_compressed_image (inputfptr, &stat) && hdutype == IMAGE_HDU &&
1140 naxis != 0 && totpix != 0 && fpvar.do_images) {
1141
1142 /* rescale a scaled integer image to reduce noise? */
1143 if (fpvar.rescale_noise != 0. && bitpix > 0 && bitpix < LONGLONG_IMG) {
1144
1145 tstatus = 0;
1146 fits_read_key(inputfptr, TDOUBLE, "BSCALE", &bscale, 0, &tstatus);
1147
1148 if (tstatus == 0 && bscale != 1.0) { /* image must be scaled */
1149
1150 if (bitpix == LONG_IMG)
1151 fp_i4stat(inputfptr, naxis, naxes, &imagestats, &stat);
1152 else
1153 fp_i2stat(inputfptr, naxis, naxes, &imagestats, &stat);
1154
1155 /* use the minimum of the MAD 2nd, 3rd, and 5th order noise estimates */
1156 noisemin = imagestats.noise3;
1157 if (imagestats.noise2 != 0. && imagestats.noise2 < noisemin) noisemin = imagestats.noise2;
1158 if (imagestats.noise5 != 0. && imagestats.noise5 < noisemin) noisemin = imagestats.noise5;
1159
1160 rescale = noisemin / fpvar.rescale_noise;
1161 if (rescale > 1.0) {
1162
1163 /* all the criteria are met, so create a temporary file that */
1164 /* contains a rescaled version of the image, in CWD */
1165
1166 /* create temporary file name */
1167 fp_tmpnam("Tmpfile3", "", tempfilename3);
1168
1169 fits_create_file(&tempfile, tempfilename3, &stat);
1170
1171 fits_get_hdu_num(inputfptr, &hdunum);
1172 if (hdunum != 1) {
1173
1174 /* the input hdu is an image extension, so create dummy primary */
1175 fits_create_img(tempfile, 8, 0, naxes, &stat);
1176 }
1177
1178 fits_copy_header(inputfptr, tempfile, &stat); /* copy the header */
1179
1180 /* rescale the data, so that it will compress more efficiently */
1181 if (bitpix == LONG_IMG)
1182 fp_i4rescale(inputfptr, naxis, naxes, rescale, tempfile, &stat);
1183 else
1184 fp_i2rescale(inputfptr, naxis, naxes, rescale, tempfile, &stat);
1185
1186 /* scale the BSCALE keyword by the inverse factor */
1187
1188 bscale = bscale * rescale;
1189 fits_update_key(tempfile, TDOUBLE, "BSCALE", &bscale, 0, &stat);
1190
1191 /* rescan the header, to reset the actual scaling parameters */
1192 fits_set_hdustruc(tempfile, &stat);
1193
1194 infptr = tempfile;
1195 rescale_flag = 1;
1196 }
1197 }
1198 }
1199
1200 if (!rescale_flag) /* just compress the input file, without rescaling */
1201 infptr = inputfptr;
1202
1203 /* compute basic statistics about the input image */
1204 if (bitpix == BYTE_IMG) {
1205 bpix = 8;
1206 strncpy(dtype, "8 ", sizeof(dtype)/sizeof(dtype[0]));
1207 fp_i2stat(infptr, naxis, naxes, &imagestats, &stat);
1208 } else if (bitpix == SHORT_IMG) {
1209 bpix = 16;
1210 strncpy(dtype, "16 ", sizeof(dtype)/sizeof(dtype[0]));
1211 fp_i2stat(infptr, naxis, naxes, &imagestats, &stat);
1212 } else if (bitpix == LONG_IMG) {
1213 bpix = 32;
1214 strncpy(dtype, "32 ", sizeof(dtype)/sizeof(dtype[0]));
1215 fp_i4stat(infptr, naxis, naxes, &imagestats, &stat);
1216 } else if (bitpix == LONGLONG_IMG) {
1217 bpix = 64;
1218 strncpy(dtype, "64 ", sizeof(dtype)/sizeof(dtype[0]));
1219 } else if (bitpix == FLOAT_IMG) {
1220 bpix = 32;
1221 strncpy(dtype, "-32", sizeof(dtype)/sizeof(dtype[0]));
1222 fp_r4stat(infptr, naxis, naxes, &imagestats, &stat);
1223 } else if (bitpix == DOUBLE_IMG) {
1224 bpix = 64;
1225 strncpy(dtype, "-64", sizeof(dtype)/sizeof(dtype[0]));
1226 fp_r4stat(infptr, naxis, naxes, &imagestats, &stat);
1227 }
1228 else dtype[0] = '\0';
1229
1230 /* use the minimum of the MAD 2nd, 3rd, and 5th order noise estimates */
1231 noisemin = imagestats.noise3;
1232 if (imagestats.noise2 != 0. && imagestats.noise2 < noisemin) noisemin = imagestats.noise2;
1233 if (imagestats.noise5 != 0. && imagestats.noise5 < noisemin) noisemin = imagestats.noise5;
1234
1235 xbits = (float) (log10(noisemin)/.301 + 1.792);
1236
1237 printf("\n File: %s\n", infits);
1238 printf(" Ext BITPIX Dimens. Nulls Min Max Mean Sigma Noise2 Noise3 Noise5 Nbits MaxR\n");
1239
1240 printf(" %3d %s", extnum, dtype);
1241 snprintf(dimen,100," (%ld", naxes[0]);
1242 len =strlen(dimen);
1243 for (ii = 1; ii < naxis; ii++) {
1244 if (len < 99)
1245 snprintf(dimen+len,100-len,",%ld", naxes[ii]);
1246 len =strlen(dimen);
1247 }
1248 if (strlen(dimen)<99)
1249 strcat(dimen, ")");
1250 printf("%-12s",dimen);
1251
1252 fits_get_hduaddr(inputfptr, &headstart, &datastart, &dataend, &stat);
1253 origdata = (float) ((dataend - datastart)/1000000.);
1254
1255 /* get elapsed and cpu times need to read the uncompressed image */
1256 fits_read_image_speed (infptr, &whole_elapse, &whole_cpu,
1257 &row_elapse, &row_cpu, &stat);
1258
1259 printf(" %5d %6.0f %6.0f %8.1f %#8.2g %#7.3g %#7.3g %#7.3g %#5.1f %#6.2f\n",
1260 imagestats.n_nulls, imagestats.minval, imagestats.maxval,
1261 imagestats.mean, imagestats.sigma,
1262 imagestats.noise2, imagestats.noise3, imagestats.noise5, xbits, bpix/xbits);
1263
1264 printf("\n Type Ratio Size (MB) Pk (Sec) UnPk Exact ElpN CPUN Elp1 CPU1\n");
1265
1266 printf(" Native %5.3f %5.3f %5.3f %5.3f\n",
1267 whole_elapse, whole_cpu, row_elapse, row_cpu);
1268
1269 if (fpvar.outfile[0]) {
1270 fprintf(outreport,
1271 " %s %d %d %ld %ld %#10.4g %d %#10.4g %#10.4g %#10.4g %#10.4g %#10.4g %#10.4g %#10.4g %#10.4g %#10.4g %#10.4g %#10.4g %#10.4g",
1272 infits, extnum, bitpix, naxes[0], naxes[1], origdata, imagestats.n_nulls, imagestats.minval,
1273 imagestats.maxval, imagestats.mean, imagestats.sigma,
1274 imagestats.noise1, imagestats.noise2, imagestats.noise3, imagestats.noise5, whole_elapse, whole_cpu, row_elapse, row_cpu);
1275 }
1276
1277 fits_set_lossy_int (outfptr, fpvar.int_to_float, &stat);
1278 if ( (bitpix > 0) && (fpvar.int_to_float != 0) ) {
1279
1280 if ( (noisemin < (fpvar.n3ratio * fpvar.quantize_level) ) ||
1281 (noisemin < fpvar.n3min)) {
1282
1283 /* image contains too little noise to quantize effectively */
1284 fits_set_lossy_int (outfptr, 0, &stat);
1285 fits_get_hdu_num(infptr, &hdunum);
1286
1287 printf(" HDU %d does not meet noise criteria to be quantized, so losslessly compressed.\n", hdunum);
1288 }
1289 }
1290
1291 /* test compression ratio and speed for each algorithm */
1292
1293 if (fpvar.quantize_level != 0) {
1294
1295 fits_set_compression_type (outfptr, RICE_1, &stat);
1296 fits_set_tile_dim (outfptr, 6, fpvar.ntile, &stat);
1297 if (fpvar.no_dither)
1298 fits_set_quantize_method(outfptr, -1, &stat);
1299 else
1300 fits_set_quantize_method(outfptr, fpvar.dither_method, &stat);
1301
1302 fits_set_quantize_level (outfptr, fpvar.quantize_level, &stat);
1303 fits_set_dither_offset(outfptr, fpvar.dither_offset, &stat);
1304 fits_set_hcomp_scale (outfptr, fpvar.scale, &stat);
1305 fits_set_hcomp_smooth (outfptr, fpvar.smooth, &stat);
1306
1307 fp_test_hdu(infptr, outfptr, outfptr2, fpvar, &stat);
1308 }
1309
1310 if (fpvar.quantize_level != 0) {
1311 \
1312 fits_set_compression_type (outfptr, HCOMPRESS_1, &stat);
1313 fits_set_tile_dim (outfptr, 6, fpvar.ntile, &stat);
1314
1315 if (fpvar.no_dither)
1316 fits_set_quantize_method(outfptr, -1, &stat);
1317 else
1318 fits_set_quantize_method(outfptr, fpvar.dither_method, &stat);
1319
1320 fits_set_quantize_level (outfptr, fpvar.quantize_level, &stat);
1321 fits_set_dither_offset(outfptr, fpvar.dither_offset, &stat);
1322 fits_set_hcomp_scale (outfptr, fpvar.scale, &stat);
1323 fits_set_hcomp_smooth (outfptr, fpvar.smooth, &stat);
1324
1325 fp_test_hdu(infptr, outfptr, outfptr2, fpvar, &stat);
1326 }
1327
1328 if (fpvar.comptype == GZIP_2) {
1329 fits_set_compression_type (outfptr, GZIP_2, &stat);
1330 } else {
1331 fits_set_compression_type (outfptr, GZIP_1, &stat);
1332 }
1333
1334 fits_set_tile_dim (outfptr, 6, fpvar.ntile, &stat);
1335
1336 if (fpvar.no_dither)
1337 fits_set_quantize_method(outfptr, -1, &stat);
1338 else
1339 fits_set_quantize_method(outfptr, fpvar.dither_method, &stat);
1340
1341 fits_set_quantize_level (outfptr, fpvar.quantize_level, &stat);
1342 fits_set_dither_offset(outfptr, fpvar.dither_offset, &stat);
1343 fits_set_hcomp_scale (outfptr, fpvar.scale, &stat);
1344 fits_set_hcomp_smooth (outfptr, fpvar.smooth, &stat);
1345
1346 fp_test_hdu(infptr, outfptr, outfptr2, fpvar, &stat);
1347
1348 /*
1349 fits_set_compression_type (outfptr, BZIP2_1, &stat);
1350 fits_set_tile_dim (outfptr, 6, fpvar.ntile, &stat);
1351 fp_test_hdu(infptr, outfptr, outfptr2, fpvar, &stat);
1352 */
1353 /*
1354 fits_set_compression_type (outfptr, PLIO_1, &stat);
1355 fits_set_tile_dim (outfptr, 6, fpvar.ntile, &stat);
1356 fp_test_hdu(infptr, outfptr, outfptr2, fpvar, &stat);
1357 */
1358 /*
1359 if (bitpix == SHORT_IMG || bitpix == LONG_IMG) {
1360 fits_set_compression_type (outfptr, NOCOMPRESS, &stat);
1361 fits_set_tile_dim (outfptr, 6, fpvar.ntile, &stat);
1362 fp_test_hdu(infptr, outfptr, outfptr2, fpvar, &stat);
1363 }
1364 */
1365 if (fpvar.outfile[0])
1366 fprintf(outreport,"\n");
1367
1368 /* delete the temporary file */
1369 if (rescale_flag) {
1370 fits_delete_file (infptr, &stat);
1371 tempfilename3[0] = '\0'; /* clear the temp filename */
1372 }
1373 } else if ( (hdutype == BINARY_TBL) && fpvar.do_tables) {
1374
1375 fits_get_num_rowsll(inputfptr, &nrows, &stat);
1376 fits_get_num_cols(inputfptr, &ncols, &stat);
1377 #if defined(_MSC_VER)
1378 /* Microsoft Visual C++ 6.0 uses '%I64d' syntax for 8-byte integers */
1379 printf("\n File: %s, HDU %d, %d cols X %I64d rows\n", infits, extnum, ncols, nrows);
1380 #elif (USE_LL_SUFFIX == 1)
1381 printf("\n File: %s, HDU %d, %d cols X %lld rows\n", infits, extnum, ncols, nrows);
1382 #else
1383 printf("\n File: %s, HDU %d, %d cols X %ld rows\n", infits, extnum, ncols, nrows);
1384 #endif
1385 fp_test_table(inputfptr, outfptr, outfptr2, fpvar, &stat);
1386
1387 } else {
1388 fits_copy_hdu (inputfptr, outfptr, 0, &stat);
1389 fits_copy_hdu (inputfptr, outfptr2, 0, &stat);
1390 }
1391
1392 fits_movrel_hdu (inputfptr, 1, NULL, &stat);
1393 extnum++;
1394 }
1395
1396
1397 if (stat == END_OF_FILE) stat = 0;
1398
1399 fits_close_file (outfptr2, &stat);
1400 fits_close_file (outfptr, &stat);
1401 fits_close_file (inputfptr, &stat);
1402
1403 if (stat) {
1404 fits_report_error (stderr, stat);
1405 }
1406 return(0);
1407 }
1408 /*--------------------------------------------------------------------------*/
fp_pack_hdu(fitsfile * infptr,fitsfile * outfptr,fpstate fpvar,int * islossless,int * status)1409 int fp_pack_hdu (fitsfile *infptr, fitsfile *outfptr, fpstate fpvar,
1410 int *islossless, int *status)
1411 {
1412 fitsfile *tempfile;
1413 long naxes[9] = {1, 1, 1, 1, 1, 1, 1, 1, 1};
1414 int stat=0, totpix=0, naxis=0, ii, hdutype, bitpix;
1415 int tstatus, hdunum;
1416 double bscale, rescale;
1417
1418 char outfits[SZ_STR], fzalgor[FLEN_VALUE];
1419 long headstart, datastart, dataend, datasize;
1420 double noisemin;
1421 /* structure to hold image statistics (defined in fpack.h) */
1422 imgstats imagestats;
1423
1424 if (*status) return(0);
1425
1426 fits_get_hdu_type (infptr, &hdutype, &stat);
1427
1428 if (hdutype == IMAGE_HDU) {
1429 fits_get_img_param (infptr, 9, &bitpix, &naxis, naxes, &stat);
1430 for (totpix=1, ii=0; ii < 9; ii++) totpix *= naxes[ii];
1431 }
1432
1433 /* check directive keyword to see if this HDU should not be compressed */
1434 tstatus = 0;
1435 if (!fits_read_key(infptr, TSTRING, "FZALGOR", fzalgor, NULL, &tstatus) ) {
1436 if (!strcmp(fzalgor, "NONE") || !strcmp(fzalgor, "none") ) {
1437 fits_copy_hdu (infptr, outfptr, 0, &stat);
1438
1439 *status = stat;
1440 return(0);
1441 }
1442 }
1443
1444 /* =============================================================== */
1445 /* This block is only for binary table compression */
1446 if (hdutype == BINARY_TBL && fpvar.do_tables) {
1447
1448 fits_get_hduaddr(infptr, &headstart, &datastart, &dataend, status);
1449 datasize = dataend - datastart;
1450
1451 if (datasize <= 2880) {
1452 /* data is less than 1 FITS block in size, so don't compress */
1453 fits_copy_hdu (infptr, outfptr, 0, &stat);
1454 } else {
1455 fits_compress_table (infptr, outfptr, &stat);
1456 }
1457
1458 *status = stat;
1459 return(0);
1460 }
1461 /* =============================================================== */
1462
1463 /* If this is not a non-null image HDU, just copy it verbatim */
1464 if (fits_is_compressed_image (infptr, &stat) || hdutype != IMAGE_HDU ||
1465 naxis == 0 || totpix == 0 || !fpvar.do_images) {
1466 fits_copy_hdu (infptr, outfptr, 0, &stat);
1467
1468 } else { /* remaining code deals only with IMAGE HDUs */
1469
1470 /* special case: rescale a scaled integer image to reduce noise? */
1471 if (fpvar.rescale_noise != 0. && bitpix > 0 && bitpix < LONGLONG_IMG) {
1472
1473 tstatus = 0;
1474 fits_read_key(infptr, TDOUBLE, "BSCALE", &bscale, 0, &tstatus);
1475 if (tstatus == 0 && bscale != 1.0) { /* image must be scaled */
1476
1477 if (bitpix == LONG_IMG)
1478 fp_i4stat(infptr, naxis, naxes, &imagestats, &stat);
1479 else
1480 fp_i2stat(infptr, naxis, naxes, &imagestats, &stat);
1481
1482 /* use the minimum of the MAD 2nd, 3rd, and 5th order noise estimates */
1483 noisemin = imagestats.noise3;
1484 if (imagestats.noise2 != 0. && imagestats.noise2 < noisemin) noisemin = imagestats.noise2;
1485 if (imagestats.noise5 != 0. && imagestats.noise5 < noisemin) noisemin = imagestats.noise5;
1486
1487 rescale = noisemin / fpvar.rescale_noise;
1488 if (rescale > 1.0) {
1489
1490 /* all the criteria are met, so create a temporary file that */
1491 /* contains a rescaled version of the image, in output directory */
1492
1493 /* create temporary file name */
1494 fits_file_name(outfptr, outfits, &stat); /* get the output file name */
1495 fp_tmpnam("Tmp3", outfits, tempfilename3);
1496
1497 fits_create_file(&tempfile, tempfilename3, &stat);
1498
1499 fits_get_hdu_num(infptr, &hdunum);
1500 if (hdunum != 1) {
1501
1502 /* the input hdu is an image extension, so create dummy primary */
1503 fits_create_img(tempfile, 8, 0, naxes, &stat);
1504 }
1505
1506 fits_copy_header(infptr, tempfile, &stat); /* copy the header */
1507
1508 /* rescale the data, so that it will compress more efficiently */
1509 if (bitpix == LONG_IMG)
1510 fp_i4rescale(infptr, naxis, naxes, rescale, tempfile, &stat);
1511 else
1512 fp_i2rescale(infptr, naxis, naxes, rescale, tempfile, &stat);
1513
1514
1515 /* scale the BSCALE keyword by the inverse factor */
1516
1517 bscale = bscale * rescale;
1518 fits_update_key(tempfile, TDOUBLE, "BSCALE", &bscale, 0, &stat);
1519
1520 /* rescan the header, to reset the actual scaling parameters */
1521 fits_set_hdustruc(tempfile, &stat);
1522
1523 fits_img_compress (tempfile, outfptr, &stat);
1524 fits_delete_file (tempfile, &stat);
1525 tempfilename3[0] = '\0'; /* clear the temp filename */
1526 *islossless = 0; /* used a lossy compression method */
1527
1528 *status = stat;
1529 return(0);
1530 }
1531 }
1532 }
1533
1534 /* if requested to do lossy compression of integer images (by */
1535 /* converting to float), then check if this HDU qualifies */
1536 if ( (bitpix > 0) && (fpvar.int_to_float != 0) ) {
1537
1538 if (bitpix >= LONG_IMG)
1539 fp_i4stat(infptr, naxis, naxes, &imagestats, &stat);
1540 else
1541 fp_i2stat(infptr, naxis, naxes, &imagestats, &stat);
1542
1543 /* rescan the image header to reset scaling values (changed by fp_iNstat) */
1544 ffrhdu(infptr, &hdutype, &stat);
1545
1546 /* use the minimum of the MAD 2nd, 3rd, and 5th order noise estimates */
1547 noisemin = imagestats.noise3;
1548 if (imagestats.noise2 != 0. && imagestats.noise2 < noisemin) noisemin = imagestats.noise2;
1549 if (imagestats.noise5 != 0. && imagestats.noise5 < noisemin) noisemin = imagestats.noise5;
1550
1551 if ( (noisemin < (fpvar.n3ratio * fpvar.quantize_level) ) ||
1552 (imagestats.noise3 < fpvar.n3min)) {
1553
1554 /* image contains too little noise to quantize effectively */
1555 fits_set_lossy_int (outfptr, 0, &stat);
1556
1557 fits_get_hdu_num(infptr, &hdunum);
1558
1559 printf(" HDU %d does not meet noise criteria to be quantized, so losslessly compressed.\n", hdunum);
1560
1561 } else {
1562 /* compressed image is not identical to original */
1563 *islossless = 0;
1564 }
1565 }
1566
1567 /* finally, do the actual image compression */
1568 fits_img_compress (infptr, outfptr, &stat);
1569
1570 if (bitpix < 0 ||
1571 (fpvar.comptype == HCOMPRESS_1 && fpvar.scale != 0.)) {
1572
1573 /* compressed image is not identical to original */
1574 *islossless = 0;
1575 }
1576 }
1577
1578 *status = stat;
1579 return(0);
1580 }
1581
1582 /*--------------------------------------------------------------------------*/
fp_unpack_hdu(fitsfile * infptr,fitsfile * outfptr,fpstate fpvar,int * status)1583 int fp_unpack_hdu (fitsfile *infptr, fitsfile *outfptr, fpstate fpvar, int *status)
1584 {
1585 UNUSED(fpvar);
1586
1587 int hdutype, lval;
1588
1589 if (*status > 0) return(0);
1590
1591 fits_get_hdu_type (infptr, &hdutype, status);
1592
1593 /* =============================================================== */
1594 /* This block is only for beta testing of binary table compression */
1595 if (hdutype == BINARY_TBL) {
1596
1597 fits_read_key(infptr, TLOGICAL, "ZTABLE", &lval, NULL, status);
1598
1599 if (*status == 0 && lval != 0) {
1600 /* uncompress the table */
1601 fits_uncompress_table (infptr, outfptr, status);
1602 } else {
1603 if (*status == KEY_NO_EXIST) /* table is not compressed */
1604 *status = 0;
1605 fits_copy_hdu (infptr, outfptr, 0, status);
1606 }
1607
1608 return(0);
1609 /* =============================================================== */
1610
1611 } else if (fits_is_compressed_image (infptr, status)) {
1612 /* uncompress the compressed image HDU */
1613 fits_img_decompress (infptr, outfptr, status);
1614 } else {
1615 /* not a compressed image HDU, so just copy it to the output */
1616 fits_copy_hdu (infptr, outfptr, 0, status);
1617 }
1618
1619 return(0);
1620 }
1621 /*--------------------------------------------------------------------------*/
fits_read_image_speed(fitsfile * infptr,float * whole_elapse,float * whole_cpu,float * row_elapse,float * row_cpu,int * status)1622 int fits_read_image_speed (fitsfile *infptr, float *whole_elapse,
1623 float *whole_cpu, float *row_elapse, float *row_cpu, int *status)
1624 {
1625 unsigned char *carray, cnull = 0;
1626 short *sarray, snull=0;
1627 int bitpix, naxis, anynull, *iarray, inull = 0;
1628 long ii, naxes[9], fpixel[9]={1,1,1,1,1,1,1,1,1}, lpixel[9]={1,1,1,1,1,1,1,1,1};
1629 long inc[9]={1,1,1,1,1,1,1,1,1} ;
1630 float *earray, enull = 0, filesize;
1631 double *darray, dnull = 0;
1632
1633 if (*status) return(*status);
1634
1635 fits_get_img_param (infptr, 9, &bitpix, &naxis, naxes, status);
1636
1637 if (naxis != 2)return(*status);
1638
1639 lpixel[0] = naxes[0];
1640 lpixel[1] = naxes[1];
1641
1642 /* filesize in MB */
1643 filesize = (float) (naxes[0] * abs(bitpix) / 8000000. * naxes[1]);
1644
1645 /* measure time required to read the raw image */
1646 fits_set_bscale(infptr, 1.0, 0.0, status);
1647 *whole_elapse = 0.;
1648 *whole_cpu = 0;
1649
1650 if (bitpix == BYTE_IMG) {
1651 carray = calloc(naxes[1]*naxes[0], sizeof(char));
1652
1653 /* remove any cached uncompressed tile
1654 (dangerous to directly modify the structure!) */
1655 /* (infptr->Fptr)->tilerow = 0; */
1656
1657 marktime(status);
1658 fits_read_subset(infptr, TBYTE, fpixel, lpixel, inc, &cnull,
1659 carray, &anynull, status);
1660
1661 /* get elapsped times */
1662 gettime(whole_elapse, whole_cpu, status);
1663
1664 /* now read the image again, row by row */
1665 if (row_elapse) {
1666
1667 /* remove any cached uncompressed tile
1668 (dangerous to directly modify the structure!) */
1669 /* (infptr->Fptr)->tilerow = 0; */
1670
1671 marktime(status);
1672 for (ii = 0; ii < naxes[1]; ii++) {
1673 fpixel[1] = ii+1;
1674 fits_read_pix(infptr, TBYTE, fpixel, naxes[0], &cnull,
1675 carray, &anynull, status);
1676 }
1677 /* get elapsped times */
1678 gettime(row_elapse, row_cpu, status);
1679 }
1680 free(carray);
1681
1682 } else if (bitpix == SHORT_IMG) {
1683 sarray = calloc(naxes[0]*naxes[1], sizeof(short));
1684
1685 marktime(status);
1686 fits_read_subset(infptr, TSHORT, fpixel, lpixel, inc, &snull,
1687 sarray, &anynull, status);
1688
1689 gettime(whole_elapse, whole_cpu, status); /* get elapsped times */
1690
1691 /* now read the image again, row by row */
1692 if (row_elapse) {
1693 marktime(status);
1694 for (ii = 0; ii < naxes[1]; ii++) {
1695
1696 fpixel[1] = ii+1;
1697 fits_read_pix(infptr, TSHORT, fpixel, naxes[0], &snull,
1698 sarray, &anynull, status);
1699 }
1700 /* get elapsped times */
1701 gettime(row_elapse, row_cpu, status);
1702 }
1703
1704 free(sarray);
1705
1706 } else if (bitpix == LONG_IMG) {
1707 iarray = calloc(naxes[0]*naxes[1], sizeof(int));
1708
1709 marktime(status);
1710
1711 fits_read_subset(infptr, TINT, fpixel, lpixel, inc, &inull,
1712 iarray, &anynull, status);
1713
1714 /* get elapsped times */
1715 gettime(whole_elapse, whole_cpu, status);
1716
1717
1718 /* now read the image again, row by row */
1719 if (row_elapse) {
1720 marktime(status);
1721 for (ii = 0; ii < naxes[1]; ii++) {
1722 fpixel[1] = ii+1;
1723 fits_read_pix(infptr, TINT, fpixel, naxes[0], &inull,
1724 iarray, &anynull, status);
1725 }
1726 /* get elapsped times */
1727 gettime(row_elapse, row_cpu, status);
1728 }
1729
1730
1731 free(iarray);
1732
1733 } else if (bitpix == FLOAT_IMG) {
1734 earray = calloc(naxes[1]*naxes[0], sizeof(float));
1735
1736 marktime(status);
1737
1738 fits_read_subset(infptr, TFLOAT, fpixel, lpixel, inc, &enull,
1739 earray, &anynull, status);
1740
1741 /* get elapsped times */
1742 gettime(whole_elapse, whole_cpu, status);
1743
1744 /* now read the image again, row by row */
1745 if (row_elapse) {
1746 marktime(status);
1747 for (ii = 0; ii < naxes[1]; ii++) {
1748 fpixel[1] = ii+1;
1749 fits_read_pix(infptr, TFLOAT, fpixel, naxes[0], &enull,
1750 earray, &anynull, status);
1751 }
1752 /* get elapsped times */
1753 gettime(row_elapse, row_cpu, status);
1754 }
1755
1756 free(earray);
1757
1758 } else if (bitpix == DOUBLE_IMG) {
1759 darray = calloc(naxes[1]*naxes[0], sizeof(double));
1760
1761 marktime(status);
1762
1763 fits_read_subset(infptr, TDOUBLE, fpixel, lpixel, inc, &dnull,
1764 darray, &anynull, status);
1765
1766 /* get elapsped times */
1767 gettime(whole_elapse, whole_cpu, status);
1768
1769 /* now read the image again, row by row */
1770 if (row_elapse) {
1771 marktime(status);
1772 for (ii = 0; ii < naxes[1]; ii++) {
1773 fpixel[1] = ii+1;
1774 fits_read_pix(infptr, TDOUBLE, fpixel, naxes[0], &dnull,
1775 darray, &anynull, status);
1776 }
1777 /* get elapsped times */
1778 gettime(row_elapse, row_cpu, status);
1779 }
1780
1781 free(darray);
1782 }
1783
1784 if (whole_elapse) *whole_elapse = *whole_elapse / filesize;
1785 if (row_elapse) *row_elapse = *row_elapse / filesize;
1786 if (whole_cpu) *whole_cpu = *whole_cpu / filesize;
1787 if (row_cpu) *row_cpu = *row_cpu / filesize;
1788
1789 return(*status);
1790 }
1791 /*--------------------------------------------------------------------------*/
fp_test_hdu(fitsfile * infptr,fitsfile * outfptr,fitsfile * outfptr2,fpstate fpvar,int * status)1792 int fp_test_hdu (fitsfile *infptr, fitsfile *outfptr, fitsfile *outfptr2,
1793 fpstate fpvar, int *status)
1794 {
1795 /* This routine is only used for performance testing of image HDUs. */
1796 /* Use fp_test_table for testing binary table HDUs. */
1797
1798 int stat = 0, hdutype, comptype;
1799 char ctype[20], lossless[4];
1800 long headstart, datastart, dataend;
1801 float origdata = 0., compressdata = 0.;
1802 float compratio = 0., packcpu = 0., unpackcpu = 0.;
1803 float elapse, whole_elapse, row_elapse, whole_cpu, row_cpu;
1804 unsigned long datasum1, datasum2, hdusum;
1805
1806 if (*status) return(0);
1807
1808 origdata = 0;
1809 compressdata = 0;
1810 compratio = 0.;
1811 lossless[0] = '\0';
1812
1813 fits_get_compression_type(outfptr, &comptype, &stat);
1814 if (comptype == RICE_1)
1815 strcpy(ctype, "RICE");
1816 else if (comptype == GZIP_1)
1817 strcpy(ctype, "GZIP1");
1818 else if (comptype == GZIP_2)
1819 strcpy(ctype, "GZIP2");/*
1820 else if (comptype == BZIP2_1)
1821 strcpy(ctype, "BZIP2");
1822 */
1823 else if (comptype == PLIO_1)
1824 strcpy(ctype, "PLIO");
1825 else if (comptype == HCOMPRESS_1)
1826 strcpy(ctype, "HCOMP");
1827 else if (comptype == NOCOMPRESS)
1828 strcpy(ctype, "NONE");
1829 else {
1830 fp_msg ("Error: unsupported image compression type ");
1831 *status = DATA_COMPRESSION_ERR;
1832 return(0);
1833 }
1834
1835 /* -------------- COMPRESS the image ------------------ */
1836
1837 marktime(&stat);
1838
1839 fits_img_compress (infptr, outfptr, &stat);
1840
1841 /* get elapsped times */
1842 gettime(&elapse, &packcpu, &stat);
1843
1844 /* get elapsed and cpu times need to read the compressed image */
1845 fits_read_image_speed (outfptr, &whole_elapse, &whole_cpu,
1846 &row_elapse, &row_cpu, &stat);
1847
1848 if (!stat) {
1849
1850 /* -------------- UNCOMPRESS the image ------------------ */
1851
1852 /* remove any cached uncompressed tile
1853 (dangerous to directly modify the structure!) */
1854 /* (outfptr->Fptr)->tilerow = 0; */
1855 marktime(&stat);
1856
1857 fits_img_decompress (outfptr, outfptr2, &stat);
1858
1859 /* get elapsped times */
1860 gettime(&elapse, &unpackcpu, &stat);
1861
1862 /* ----------------------------------------------------- */
1863
1864
1865 /* get sizes of original and compressed images */
1866
1867 fits_get_hduaddr(infptr, &headstart, &datastart, &dataend, &stat);
1868 origdata = (float) ((dataend - datastart)/1000000.);
1869
1870 fits_get_hduaddr(outfptr, &headstart, &datastart, &dataend, &stat);
1871 compressdata = (float) ((dataend - datastart)/1000000.);
1872
1873 if (compressdata != 0)
1874 compratio = (float) origdata / (float) compressdata;
1875
1876 /* is this uncompressed image identical to the original? */
1877
1878 fits_get_chksum(infptr, &datasum1, &hdusum, &stat);
1879 fits_get_chksum(outfptr2, &datasum2, &hdusum, &stat);
1880
1881 if ( datasum1 == datasum2) {
1882 strcpy(lossless, "Yes");
1883 } else {
1884 strcpy(lossless, "No");
1885 }
1886
1887 printf(" %-5s %6.2f %7.2f ->%7.2f %7.2f %7.2f %s %5.3f %5.3f %5.3f %5.3f\n",
1888 ctype, compratio, origdata, compressdata,
1889 packcpu, unpackcpu, lossless, whole_elapse, whole_cpu,
1890 row_elapse, row_cpu);
1891
1892
1893 if (fpvar.outfile[0]) {
1894 fprintf(outreport," %6.3f %5.2f %5.2f %s %7.3f %7.3f %7.3f %7.3f",
1895 compratio, packcpu, unpackcpu, lossless, whole_elapse, whole_cpu,
1896 row_elapse, row_cpu);
1897 }
1898
1899 /* delete the output HDUs to concerve disk space */
1900
1901 fits_delete_hdu(outfptr, &hdutype, &stat);
1902 fits_delete_hdu(outfptr2, &hdutype, &stat);
1903
1904 } else {
1905 printf(" %-5s (unable to compress image)\n", ctype);
1906 }
1907
1908 /* try to recover from any compression errors */
1909 if (stat == DATA_COMPRESSION_ERR) stat = 0;
1910
1911 *status = stat;
1912 return(0);
1913 }
1914 /*--------------------------------------------------------------------------*/
fp_test_table(fitsfile * infptr,fitsfile * outfptr,fitsfile * outfptr2,fpstate fpvar,int * status)1915 int fp_test_table (fitsfile *infptr, fitsfile *outfptr, fitsfile *outfptr2,
1916 fpstate fpvar, int *status)
1917 {
1918 /* this routine is for performance testing of the table compression methods */
1919 UNUSED(fpvar);
1920 UNUSED(outfptr2);
1921
1922 int stat = 0, hdutype, tstatus = 0;
1923 char fzalgor[FLEN_VALUE];
1924 LONGLONG headstart, datastart, dataend;
1925 float elapse, cpu;
1926
1927 if (*status) return(0);
1928
1929 /* check directive keyword to see if this HDU should not be compressed */
1930 if (!fits_read_key(infptr, TSTRING, "FZALGOR", fzalgor, NULL, &tstatus) ) {
1931 if (!strcmp(fzalgor, "NONE") || !strcmp(fzalgor, "none")) {
1932 return(0);
1933 }
1934 }
1935
1936 fits_get_hduaddrll(infptr, &headstart, &datastart, &dataend, status);
1937
1938 /* can't compress small tables with less than 2880 bytes of data */
1939 if (dataend - datastart <= 2880) {
1940 return(0);
1941 }
1942
1943 marktime(&stat);
1944 stat= -999; /* set special flag value */
1945 fits_compress_table (infptr, outfptr, &stat);
1946
1947 /* get elapsped times */
1948 gettime(&elapse, &cpu, &stat);
1949
1950 fits_delete_hdu(outfptr, &hdutype, &stat);
1951
1952 printf("\nElapsed time = %f, cpu = %f\n", elapse, cpu);
1953
1954 fits_report_error (stderr, stat);
1955
1956 return(0);
1957 }
1958 /*--------------------------------------------------------------------------*/
marktime(int * status)1959 int marktime(int *status)
1960 {
1961 #if defined(unix) || defined(__unix__) || defined(__unix)
1962 struct timeval tv;
1963 /* struct timezone tz; */
1964
1965 /* gettimeofday (&tv, &tz); */
1966 gettimeofday (&tv, NULL);
1967
1968 startsec = tv.tv_sec;
1969 startmilli = tv.tv_usec/1000;
1970
1971 scpu = clock();
1972 #else
1973 /* don't support high timing precision on Windows machines */
1974 startsec = 0;
1975 startmilli = 0;
1976
1977 scpu = clock();
1978 #endif
1979 return( *status );
1980 }
1981 /*--------------------------------------------------------------------------*/
gettime(float * elapse,float * elapscpu,int * status)1982 int gettime(float *elapse, float *elapscpu, int *status)
1983 {
1984 #if defined(unix) || defined(__unix__) || defined(__unix)
1985 struct timeval tv;
1986 /* struct timezone tz; */
1987 int stopmilli;
1988 long stopsec;
1989
1990 /* gettimeofday (&tv, &tz); */
1991 gettimeofday (&tv, NULL);
1992 ecpu = clock();
1993
1994 stopmilli = tv.tv_usec/1000;
1995 stopsec = tv.tv_sec;
1996
1997 *elapse = (stopsec - startsec) + (stopmilli - startmilli)/1000.;
1998 *elapscpu = (ecpu - scpu) * 1.0 / CLOCKTICKS;
1999 /*
2000 printf(" (start: %ld + %d), stop: (%ld + %d) elapse: %f\n ",
2001 startsec,startmilli,stopsec, stopmilli, *elapse);
2002 */
2003 #else
2004 /* set the elapsed time the same as the CPU time on Windows machines */
2005 *elapscpu = (float) ((ecpu - scpu) * 1.0 / CLOCKTICKS);
2006 *elapse = *elapscpu;
2007 #endif
2008 return( *status );
2009 }
2010 /*--------------------------------------------------------------------------*/
fp_i2stat(fitsfile * infptr,int naxis,long * naxes,imgstats * imagestats,int * status)2011 int fp_i2stat(fitsfile *infptr, int naxis, long *naxes, imgstats *imagestats, int *status)
2012 {
2013 /*
2014 read the central XSAMPLE by YSAMPLE region of pixels in the int*2 image,
2015 and then compute basic statistics: min, max, mean, sigma, mean diff, etc.
2016 */
2017
2018 long fpixel[9] = {1,1,1,1,1,1,1,1,1};
2019 long lpixel[9] = {1,1,1,1,1,1,1,1,1};
2020 long inc[9] = {1,1,1,1,1,1,1,1,1};
2021 long i1, i2, npix, ngood, nx, ny;
2022 short *intarray, minvalue, maxvalue, nullvalue;
2023 int anynul, tstatus, checknull = 1;
2024 double mean, sigma, noise1, noise2, noise3, noise5;
2025
2026 /* select the middle XSAMPLE by YSAMPLE area of the image */
2027 i1 = naxes[0]/2 - (XSAMPLE/2 - 1);
2028 i2 = naxes[0]/2 + (XSAMPLE/2);
2029 if (i1 < 1) i1 = 1;
2030 if (i2 > naxes[0]) i2 = naxes[0];
2031 fpixel[0] = i1;
2032 lpixel[0] = i2;
2033 nx = i2 - i1 +1;
2034
2035 if (naxis > 1) {
2036 i1 = naxes[1]/2 - (YSAMPLE/2 - 1);
2037 i2 = naxes[1]/2 + (YSAMPLE/2);
2038 if (i1 < 1) i1 = 1;
2039 if (i2 > naxes[1]) i2 = naxes[1];
2040 fpixel[1] = i1;
2041 lpixel[1] = i2;
2042 }
2043 ny = i2 - i1 +1;
2044
2045 npix = nx * ny;
2046
2047 /* if there are higher dimensions, read the middle plane of the cube */
2048 if (naxis > 2) {
2049 fpixel[2] = naxes[2]/2 + 1;
2050 lpixel[2] = naxes[2]/2 + 1;
2051 }
2052
2053 intarray = calloc(npix, sizeof(short));
2054 if (!intarray) {
2055 *status = MEMORY_ALLOCATION;
2056 return(*status);
2057 }
2058
2059 /* turn off any scaling of the integer pixel values */
2060 fits_set_bscale(infptr, 1.0, 0.0, status);
2061
2062 fits_read_subset_sht(infptr, 0, naxis, naxes, fpixel, lpixel, inc,
2063 0, intarray, &anynul, status);
2064
2065 /* read the null value keyword (BLANK) if present */
2066 tstatus = 0;
2067 fits_read_key(infptr, TSHORT, "BLANK", &nullvalue, 0, &tstatus);
2068 if (tstatus) {
2069 nullvalue = 0;
2070 checknull = 0;
2071 }
2072
2073 /* compute statistics of the image */
2074
2075 fits_img_stats_short(intarray, nx, ny, checknull, nullvalue,
2076 &ngood, &minvalue, &maxvalue, &mean, &sigma, &noise1, &noise2, &noise3, &noise5, status);
2077
2078 imagestats->n_nulls = npix - ngood;
2079 imagestats->minval = minvalue;
2080 imagestats->maxval = maxvalue;
2081 imagestats->mean = mean;
2082 imagestats->sigma = sigma;
2083 imagestats->noise1 = noise1;
2084 imagestats->noise2 = noise2;
2085 imagestats->noise3 = noise3;
2086 imagestats->noise5 = noise5;
2087
2088 free(intarray);
2089 return(*status);
2090 }
2091 /*--------------------------------------------------------------------------*/
fp_i4stat(fitsfile * infptr,int naxis,long * naxes,imgstats * imagestats,int * status)2092 int fp_i4stat(fitsfile *infptr, int naxis, long *naxes, imgstats *imagestats, int *status)
2093 {
2094 /*
2095 read the central XSAMPLE by YSAMPLE region of pixels in the int*2 image,
2096 and then compute basic statistics: min, max, mean, sigma, mean diff, etc.
2097 */
2098
2099 long fpixel[9] = {1,1,1,1,1,1,1,1,1};
2100 long lpixel[9] = {1,1,1,1,1,1,1,1,1};
2101 long inc[9] = {1,1,1,1,1,1,1,1,1};
2102 long i1, i2, npix, ngood, nx, ny;
2103 int *intarray, minvalue, maxvalue, nullvalue;
2104 int anynul, tstatus, checknull = 1;
2105 double mean, sigma, noise1, noise2, noise3, noise5;
2106
2107 /* select the middle XSAMPLE by YSAMPLE area of the image */
2108 i1 = naxes[0]/2 - (XSAMPLE/2 - 1);
2109 i2 = naxes[0]/2 + (XSAMPLE/2);
2110 if (i1 < 1) i1 = 1;
2111 if (i2 > naxes[0]) i2 = naxes[0];
2112 fpixel[0] = i1;
2113 lpixel[0] = i2;
2114 nx = i2 - i1 +1;
2115
2116 if (naxis > 1) {
2117 i1 = naxes[1]/2 - (YSAMPLE/2 - 1);
2118 i2 = naxes[1]/2 + (YSAMPLE/2);
2119 if (i1 < 1) i1 = 1;
2120 if (i2 > naxes[1]) i2 = naxes[1];
2121 fpixel[1] = i1;
2122 lpixel[1] = i2;
2123 }
2124 ny = i2 - i1 +1;
2125
2126 npix = nx * ny;
2127
2128 /* if there are higher dimensions, read the middle plane of the cube */
2129 if (naxis > 2) {
2130 fpixel[2] = naxes[2]/2 + 1;
2131 lpixel[2] = naxes[2]/2 + 1;
2132 }
2133
2134 intarray = calloc(npix, sizeof(int));
2135 if (!intarray) {
2136 *status = MEMORY_ALLOCATION;
2137 return(*status);
2138 }
2139
2140 /* turn off any scaling of the integer pixel values */
2141 fits_set_bscale(infptr, 1.0, 0.0, status);
2142
2143 fits_read_subset_int(infptr, 0, naxis, naxes, fpixel, lpixel, inc,
2144 0, intarray, &anynul, status);
2145
2146 /* read the null value keyword (BLANK) if present */
2147 tstatus = 0;
2148 fits_read_key(infptr, TINT, "BLANK", &nullvalue, 0, &tstatus);
2149 if (tstatus) {
2150 nullvalue = 0;
2151 checknull = 0;
2152 }
2153
2154 /* compute statistics of the image */
2155
2156 fits_img_stats_int(intarray, nx, ny, checknull, nullvalue,
2157 &ngood, &minvalue, &maxvalue, &mean, &sigma, &noise1, &noise2, &noise3, &noise5, status);
2158
2159 imagestats->n_nulls = npix - ngood;
2160 imagestats->minval = minvalue;
2161 imagestats->maxval = maxvalue;
2162 imagestats->mean = mean;
2163 imagestats->sigma = sigma;
2164 imagestats->noise1 = noise1;
2165 imagestats->noise2 = noise2;
2166 imagestats->noise3 = noise3;
2167 imagestats->noise5 = noise5;
2168
2169 free(intarray);
2170 return(*status);
2171 }
2172 /*--------------------------------------------------------------------------*/
fp_r4stat(fitsfile * infptr,int naxis,long * naxes,imgstats * imagestats,int * status)2173 int fp_r4stat(fitsfile *infptr, int naxis, long *naxes, imgstats *imagestats, int *status)
2174 {
2175 /*
2176 read the central XSAMPLE by YSAMPLE region of pixels in the int*2 image,
2177 and then compute basic statistics: min, max, mean, sigma, mean diff, etc.
2178 */
2179
2180 long fpixel[9] = {1,1,1,1,1,1,1,1,1};
2181 long lpixel[9] = {1,1,1,1,1,1,1,1,1};
2182 long inc[9] = {1,1,1,1,1,1,1,1,1};
2183 long i1, i2, npix, ngood, nx, ny;
2184 float *array, minvalue, maxvalue, nullvalue = FLOATNULLVALUE;
2185 int anynul,checknull = 1;
2186 double mean, sigma, noise1, noise2, noise3, noise5;
2187
2188 /* select the middle XSAMPLE by YSAMPLE area of the image */
2189 i1 = naxes[0]/2 - (XSAMPLE/2 - 1);
2190 i2 = naxes[0]/2 + (XSAMPLE/2);
2191 if (i1 < 1) i1 = 1;
2192 if (i2 > naxes[0]) i2 = naxes[0];
2193 fpixel[0] = i1;
2194 lpixel[0] = i2;
2195 nx = i2 - i1 +1;
2196
2197 if (naxis > 1) {
2198 i1 = naxes[1]/2 - (YSAMPLE/2 - 1);
2199 i2 = naxes[1]/2 + (YSAMPLE/2);
2200 if (i1 < 1) i1 = 1;
2201 if (i2 > naxes[1]) i2 = naxes[1];
2202 fpixel[1] = i1;
2203 lpixel[1] = i2;
2204 }
2205 ny = i2 - i1 +1;
2206
2207 npix = nx * ny;
2208
2209 /* if there are higher dimensions, read the middle plane of the cube */
2210 if (naxis > 2) {
2211 fpixel[2] = naxes[2]/2 + 1;
2212 lpixel[2] = naxes[2]/2 + 1;
2213 }
2214
2215 array = calloc(npix, sizeof(float));
2216 if (!array) {
2217 *status = MEMORY_ALLOCATION;
2218 return(*status);
2219 }
2220
2221 fits_read_subset_flt(infptr, 0, naxis, naxes, fpixel, lpixel, inc,
2222 nullvalue, array, &anynul, status);
2223
2224 /* are there any null values in the array? */
2225 if (!anynul) {
2226 nullvalue = 0.;
2227 checknull = 0;
2228 }
2229
2230 /* compute statistics of the image */
2231
2232 fits_img_stats_float(array, nx, ny, checknull, nullvalue,
2233 &ngood, &minvalue, &maxvalue, &mean, &sigma, &noise1, &noise2, &noise3, &noise5, status);
2234
2235 imagestats->n_nulls = npix - ngood;
2236 imagestats->minval = minvalue;
2237 imagestats->maxval = maxvalue;
2238 imagestats->mean = mean;
2239 imagestats->sigma = sigma;
2240 imagestats->noise1 = noise1;
2241 imagestats->noise2 = noise2;
2242 imagestats->noise3 = noise3;
2243 imagestats->noise5 = noise5;
2244
2245 free(array);
2246 return(*status);
2247 }
2248 /*--------------------------------------------------------------------------*/
fp_i2rescale(fitsfile * infptr,int naxis,long * naxes,double rescale,fitsfile * outfptr,int * status)2249 int fp_i2rescale(fitsfile *infptr, int naxis, long *naxes, double rescale,
2250 fitsfile *outfptr, int *status)
2251 {
2252 /*
2253 divide the integer pixel values in the input file by rescale,
2254 and write back out to the output file..
2255 */
2256
2257 long ii, jj, nelem = 1, nx, ny;
2258 short *intarray, nullvalue;
2259 int anynul, tstatus, checknull = 1;
2260
2261 nx = naxes[0];
2262 ny = 1;
2263
2264 for (ii = 1; ii < naxis; ii++) {
2265 ny = ny * naxes[ii];
2266 }
2267
2268 intarray = calloc(nx, sizeof(short));
2269 if (!intarray) {
2270 *status = MEMORY_ALLOCATION;
2271 return(*status);
2272 }
2273
2274 /* read the null value keyword (BLANK) if present */
2275 tstatus = 0;
2276 fits_read_key(infptr, TSHORT, "BLANK", &nullvalue, 0, &tstatus);
2277 if (tstatus) {
2278 checknull = 0;
2279 }
2280
2281 /* turn off any scaling of the integer pixel values */
2282 fits_set_bscale(infptr, 1.0, 0.0, status);
2283 fits_set_bscale(outfptr, 1.0, 0.0, status);
2284
2285 for (ii = 0; ii < ny; ii++) {
2286
2287 fits_read_img_sht(infptr, 1, nelem, nx,
2288 0, intarray, &anynul, status);
2289
2290 if (checknull) {
2291 for (jj = 0; jj < nx; jj++) {
2292 if (intarray[jj] != nullvalue)
2293 intarray[jj] = NSHRT( (intarray[jj] / rescale) );
2294 }
2295 } else {
2296 for (jj = 0; jj < nx; jj++)
2297 intarray[jj] = NSHRT( (intarray[jj] / rescale) );
2298 }
2299
2300 fits_write_img_sht(outfptr, 1, nelem, nx, intarray, status);
2301
2302 nelem += nx;
2303 }
2304
2305 free(intarray);
2306 return(*status);
2307 }
2308 /*--------------------------------------------------------------------------*/
fp_i4rescale(fitsfile * infptr,int naxis,long * naxes,double rescale,fitsfile * outfptr,int * status)2309 int fp_i4rescale(fitsfile *infptr, int naxis, long *naxes, double rescale,
2310 fitsfile *outfptr, int *status)
2311 {
2312 /*
2313 divide the integer pixel values in the input file by rescale,
2314 and write back out to the output file..
2315 */
2316
2317 long ii, jj, nelem = 1, nx, ny;
2318 int *intarray, nullvalue;
2319 int anynul, tstatus, checknull = 1;
2320
2321 nx = naxes[0];
2322 ny = 1;
2323
2324 for (ii = 1; ii < naxis; ii++) {
2325 ny = ny * naxes[ii];
2326 }
2327
2328 intarray = calloc(nx, sizeof(int));
2329 if (!intarray) {
2330 *status = MEMORY_ALLOCATION;
2331 return(*status);
2332 }
2333
2334 /* read the null value keyword (BLANK) if present */
2335 tstatus = 0;
2336 fits_read_key(infptr, TINT, "BLANK", &nullvalue, 0, &tstatus);
2337 if (tstatus) {
2338 checknull = 0;
2339 }
2340
2341 /* turn off any scaling of the integer pixel values */
2342 fits_set_bscale(infptr, 1.0, 0.0, status);
2343 fits_set_bscale(outfptr, 1.0, 0.0, status);
2344
2345 for (ii = 0; ii < ny; ii++) {
2346
2347 fits_read_img_int(infptr, 1, nelem, nx,
2348 0, intarray, &anynul, status);
2349
2350 if (checknull) {
2351 for (jj = 0; jj < nx; jj++) {
2352 if (intarray[jj] != nullvalue)
2353 intarray[jj] = NINT( (intarray[jj] / rescale) );
2354 }
2355 } else {
2356 for (jj = 0; jj < nx; jj++)
2357 intarray[jj] = NINT( (intarray[jj] / rescale) );
2358 }
2359
2360 fits_write_img_int(outfptr, 1, nelem, nx, intarray, status);
2361
2362 nelem += nx;
2363 }
2364
2365 free(intarray);
2366 return(*status);
2367 }
2368 /* ========================================================================
2369 * Signal and error handler.
2370 */
abort_fpack(int sig)2371 void abort_fpack(int sig)
2372 {
2373 /* clean up by deleting temporary files */
2374 UNUSED(sig);
2375
2376 if (tempfilename[0]) {
2377 remove(tempfilename);
2378 }
2379 if (tempfilename2[0]) {
2380 remove(tempfilename2);
2381 }
2382 if (tempfilename3[0]) {
2383 remove(tempfilename3);
2384 }
2385 exit(-1);
2386 }
2387