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