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