File indexing completed on 2024-03-24 03:46:42
0001 /* FPACK utility routines 0002 SPDX-FileCopyrightText: William D. Pence <https://heasarc.gsfc.nasa.gov/fitsio/> 0003 SPDX-FileCopyrightText: R. Seaman 0004 0005 SPDX-License-Identifier: LicenseRef-NASA-FV-License-Agreement 0006 */ 0007 0008 /* The following unused functions have been deleted to suppress compiler warnings: 0009 * int fp_pack_data_to_fits (const char *inputBuffer, size_t inputBufferSize, fitsfile **outfits, fpstate fpvar, int *islossless) 0010 * int fp_pack_data_to_data (const char *inputBuffer, size_t inputBufferSize, unsigned char **outputBuffer, size_t *outputBufferSize, 0011 * fpstate fpvar, int *islossless) 0012 * int fp_unpack_data_to_fits (const char *inputBuffer, size_t inputBufferSize, fitsfile **outfits, fpstate fpvar) 0013 * 0014 * If any of these are required in the future, please refer back to before commit 56a544be 1st June 2022 0015 */ 0016 0017 #include <time.h> 0018 #include <float.h> 0019 #include <signal.h> 0020 #include <ctype.h> 0021 0022 /* #include "bzlib.h" only for experimental purposes */ 0023 0024 #if defined(unix) || defined(__unix__) || defined(__unix) 0025 #include <sys/time.h> 0026 #endif 0027 0028 #include <math.h> 0029 #include <fitsio.h> 0030 #include <fitsio2.h> 0031 #include "fpack.h" 0032 0033 /* these filename buffer are used to delete temporary files */ 0034 /* in case the program is aborted */ 0035 char tempfilename[SZ_STR]; 0036 char tempfilename2[SZ_STR]; 0037 char tempfilename3[SZ_STR]; 0038 0039 /* nearest integer function */ 0040 # define NINT(x) ((x >= 0.) ? (int) (x + 0.5) : (int) (x - 0.5)) 0041 # define NSHRT(x) ((x >= 0.) ? (short) (x + 0.5) : (short) (x - 0.5)) 0042 0043 /* define variables for measuring elapsed time */ 0044 clock_t scpu, ecpu; 0045 long startsec; /* start of elapsed time interval */ 0046 int startmilli; /* start of elapsed time interval */ 0047 0048 /* CLOCKS_PER_SEC should be defined by most compilers */ 0049 #if defined(CLOCKS_PER_SEC) 0050 #define CLOCKTICKS CLOCKS_PER_SEC 0051 #else 0052 /* on SUN OS machine, CLOCKS_PER_SEC is not defined, so set its value */ 0053 #define CLOCKTICKS 1000000 0054 #endif 0055 0056 FILE *outreport; 0057 0058 /* dimension of central image area to be sampled for test statistics */ 0059 int XSAMPLE = 4100; 0060 int YSAMPLE = 4100; 0061 0062 #define UNUSED(x) (void)(x) 0063 #define fp_tmpnam(suffix, rootname, tmpnam) _fp_tmpnam((char *)suffix, (char *)rootname, (char *)tmpnam) 0064 0065 /*--------------------------------------------------------------------------*/ 0066 int fp_noop (void) 0067 { 0068 fp_msg ("Input and output files are unchanged.\n"); 0069 return(0); 0070 } 0071 /*--------------------------------------------------------------------------*/ 0072 void fp_abort_output (fitsfile *infptr, fitsfile *outfptr, int stat) 0073 { 0074 int status = 0, hdunum; 0075 char msg[SZ_STR]; 0076 0077 if (infptr) 0078 { 0079 fits_file_name(infptr, tempfilename, &status); 0080 fits_get_hdu_num(infptr, &hdunum); 0081 0082 fits_close_file (infptr, &status); 0083 0084 snprintf(msg, SZ_STR,"Error processing file: %s\n", tempfilename); 0085 fp_msg (msg); 0086 snprintf(msg, SZ_STR," in HDU number %d\n", hdunum); 0087 fp_msg (msg); 0088 } 0089 else 0090 { 0091 snprintf(msg, SZ_STR,"Error: Unable to process input file\n"); 0092 fp_msg(msg); 0093 } 0094 fits_report_error (stderr, stat); 0095 0096 if (outfptr) { 0097 fits_delete_file(outfptr, &status); 0098 fp_msg ("Input file is unchanged.\n"); 0099 } 0100 } 0101 /*--------------------------------------------------------------------------*/ 0102 int fp_version (void) 0103 { 0104 float version; 0105 char cfitsioversion[40]; 0106 0107 fp_msg (FPACK_VERSION); 0108 fits_get_version(&version); 0109 snprintf(cfitsioversion, 40," CFITSIO version %5.3f", version); 0110 fp_msg(cfitsioversion); 0111 fp_msg ("\n"); 0112 return(0); 0113 } 0114 /*--------------------------------------------------------------------------*/ 0115 int fp_access (char *filename) 0116 { 0117 /* test if a file exists */ 0118 0119 FILE *diskfile; 0120 0121 diskfile = fopen(filename, "r"); 0122 0123 if (diskfile) { 0124 fclose(diskfile); 0125 return(0); 0126 } else { 0127 return(-1); 0128 } 0129 } 0130 /*--------------------------------------------------------------------------*/ 0131 int _fp_tmpnam(char *suffix, char *rootname, char *tmpnam) 0132 { 0133 /* create temporary file name */ 0134 0135 int maxtry = 30, ii; 0136 0137 if (strlen(suffix) + strlen(rootname) > SZ_STR-5) { 0138 fp_msg ("Error: filename is too long to create temporary file\n"); exit (-1); 0139 } 0140 0141 strcpy (tmpnam, rootname); /* start with rootname */ 0142 strcat(tmpnam, suffix); /* append the suffix */ 0143 0144 maxtry = SZ_STR - strlen(tmpnam) - 1; 0145 0146 for (ii = 0; ii < maxtry; ii++) { 0147 if (fp_access(tmpnam)) break; /* good, the file does not exist */ 0148 if (strlen(tmpnam) > SZ_STR-2) 0149 { 0150 fp_msg ("\nCould not create temporary file name:\n"); 0151 fp_msg (tmpnam); 0152 fp_msg ("\n"); 0153 exit (-1); 0154 } 0155 strcat(tmpnam, "x"); /* append an x to the name, and try again */ 0156 } 0157 0158 if (ii == maxtry) { 0159 fp_msg ("\nCould not create temporary file name:\n"); 0160 fp_msg (tmpnam); 0161 fp_msg ("\n"); 0162 exit (-1); 0163 } 0164 0165 return(0); 0166 } 0167 /*--------------------------------------------------------------------------*/ 0168 int fp_init (fpstate *fpptr) 0169 { 0170 int ii; 0171 0172 fpptr->comptype = RICE_1; 0173 fpptr->quantize_level = DEF_QLEVEL; 0174 fpptr->no_dither = 0; 0175 fpptr->dither_method = 1; 0176 fpptr->dither_offset = 0; 0177 fpptr->int_to_float = 0; 0178 0179 /* thresholds when using the -i2f flag */ 0180 fpptr->n3ratio = 2.0; /* minimum ratio of image noise sigma / q */ 0181 fpptr->n3min = 6.; /* minimum noise sigma. */ 0182 0183 fpptr->scale = DEF_HCOMP_SCALE; 0184 fpptr->smooth = DEF_HCOMP_SMOOTH; 0185 fpptr->rescale_noise = DEF_RESCALE_NOISE; 0186 fpptr->ntile[0] = (long) -1; /* -1 means extent of axis */ 0187 0188 for (ii=1; ii < MAX_COMPRESS_DIM; ii++) 0189 fpptr->ntile[ii] = (long) 1; 0190 0191 fpptr->to_stdout = 0; 0192 fpptr->listonly = 0; 0193 fpptr->clobber = 0; 0194 fpptr->delete_input = 0; 0195 fpptr->do_not_prompt = 0; 0196 fpptr->do_checksums = 1; 0197 fpptr->do_gzip_file = 0; 0198 fpptr->do_tables = 0; /* this is intended for testing purposes */ 0199 fpptr->do_images = 1; /* can be turned off with -tableonly switch */ 0200 fpptr->test_all = 0; 0201 fpptr->verbose = 0; 0202 0203 fpptr->prefix[0] = 0; 0204 fpptr->extname[0] = 0; 0205 fpptr->delete_suffix = 0; 0206 fpptr->outfile[0] = 0; 0207 0208 fpptr->firstfile = 1; 0209 0210 /* magic number for initialization check, boolean for preflight 0211 */ 0212 fpptr->initialized = FP_INIT_MAGIC; 0213 fpptr->preflight_checked = 0; 0214 return(0); 0215 } 0216 /*--------------------------------------------------------------------------*/ 0217 int fp_list (int argc, char *argv[], fpstate fpvar) 0218 { 0219 fitsfile *infptr; 0220 char infits[SZ_STR], msg[SZ_STR]; 0221 int hdunum, iarg, stat=0; 0222 LONGLONG sizell; 0223 0224 if (fpvar.initialized != FP_INIT_MAGIC) { 0225 fp_msg ("Error: internal initialization error\n"); exit (-1); 0226 } 0227 0228 for (iarg=fpvar.firstfile; iarg < argc; iarg++) { 0229 strncpy (infits, argv[iarg], SZ_STR-1); 0230 infits[SZ_STR-1]=0; 0231 0232 if (strchr (infits, '[') || strchr (infits, ']')) { 0233 fp_msg ("Error: section/extension notation not supported: "); 0234 fp_msg (infits); fp_msg ("\n"); exit (-1); 0235 } 0236 0237 if (fp_access (infits) != 0) { 0238 fp_msg ("Error: can't find or read input file "); fp_msg (infits); 0239 fp_msg ("\n"); fp_noop (); exit (-1); 0240 } 0241 0242 fits_open_file (&infptr, infits, READONLY, &stat); 0243 if (stat) { fits_report_error (stderr, stat); exit (stat); } 0244 0245 /* move to the end of file, to get the total size in bytes */ 0246 fits_get_num_hdus (infptr, &hdunum, &stat); 0247 fits_movabs_hdu (infptr, hdunum, NULL, &stat); 0248 fits_get_hduaddrll(infptr, NULL, NULL, &sizell, &stat); 0249 0250 0251 if (stat) { 0252 fp_abort_output(infptr, NULL, stat); 0253 } 0254 0255 snprintf (msg, SZ_STR,"# %s (", infits); fp_msg (msg); 0256 0257 #if defined(_MSC_VER) 0258 /* Microsoft Visual C++ 6.0 uses '%I64d' syntax for 8-byte integers */ 0259 snprintf(msg, SZ_STR,"%I64d bytes)\n", sizell); fp_msg (msg); 0260 #elif (USE_LL_SUFFIX == 1) 0261 snprintf(msg, SZ_STR,"%lld bytes)\n", sizell); fp_msg (msg); 0262 #else 0263 snprintf(msg, SZ_STR,"%ld bytes)\n", sizell); fp_msg (msg); 0264 #endif 0265 fp_info_hdu (infptr); 0266 0267 fits_close_file (infptr, &stat); 0268 if (stat) { fits_report_error (stderr, stat); exit (stat); } 0269 } 0270 return(0); 0271 } 0272 /*--------------------------------------------------------------------------*/ 0273 int fp_info_hdu (fitsfile *infptr) 0274 { 0275 long naxes[9] = {1, 1, 1, 1, 1, 1, 1, 1, 1}; 0276 char msg[SZ_STR], val[SZ_CARD], com[SZ_CARD]; 0277 int naxis=0, hdutype, bitpix, hdupos, stat=0, ii; 0278 unsigned long datasum, hdusum; 0279 0280 fits_movabs_hdu (infptr, 1, NULL, &stat); 0281 if (stat) { 0282 fp_abort_output(infptr, NULL, stat); 0283 } 0284 0285 for (hdupos=1; ! stat; hdupos++) { 0286 fits_get_hdu_type (infptr, &hdutype, &stat); 0287 if (stat) { 0288 fp_abort_output(infptr, NULL, stat); 0289 } 0290 0291 /* fits_get_hdu_type calls unknown extensions "IMAGE_HDU" 0292 * so consult XTENSION keyword itself 0293 */ 0294 fits_read_keyword (infptr, "XTENSION", val, com, &stat); 0295 if (stat == KEY_NO_EXIST) { 0296 /* in primary HDU which by definition is an "image" */ 0297 stat=0; /* clear for later error handling */ 0298 0299 } else if (stat) { 0300 fp_abort_output(infptr, NULL, stat); 0301 0302 } else if (hdutype == IMAGE_HDU) { 0303 /* that is, if XTENSION != "IMAGE" AND != "BINTABLE" */ 0304 if (strncmp (val+1, "IMAGE", 5) && 0305 strncmp (val+1, "BINTABLE", 5)) { 0306 0307 /* assign something other than any of these */ 0308 hdutype = IMAGE_HDU + ASCII_TBL + BINARY_TBL; 0309 } 0310 } 0311 0312 fits_get_chksum(infptr, &datasum, &hdusum, &stat); 0313 0314 if (hdutype == IMAGE_HDU) { 0315 snprintf (msg, SZ_STR," %d IMAGE", hdupos); fp_msg (msg); 0316 snprintf (msg, SZ_STR," SUMS=%lu/%lu", (unsigned long) (~((int) hdusum)), datasum); fp_msg (msg); 0317 0318 fits_get_img_param (infptr, 9, &bitpix, &naxis, naxes, &stat); 0319 0320 snprintf (msg, SZ_STR," BITPIX=%d", bitpix); fp_msg (msg); 0321 0322 if (naxis == 0) { 0323 snprintf (msg, SZ_STR," [no_pixels]"); fp_msg (msg); 0324 } else if (naxis == 1) { 0325 snprintf (msg, SZ_STR," [%ld]", naxes[1]); fp_msg (msg); 0326 } else { 0327 snprintf (msg, SZ_STR," [%ld", naxes[0]); fp_msg (msg); 0328 for (ii=1; ii < naxis; ii++) { 0329 snprintf (msg, SZ_STR,"x%ld", naxes[ii]); fp_msg (msg); 0330 } 0331 fp_msg ("]"); 0332 } 0333 0334 if (fits_is_compressed_image (infptr, &stat)) { 0335 fits_read_keyword (infptr, "ZCMPTYPE", val, com, &stat); 0336 0337 /* allow for quote in keyword value */ 0338 if (! strncmp (val+1, "RICE_1", 6)) 0339 fp_msg (" tiled_rice\n"); 0340 else if (! strncmp (val+1, "GZIP_1", 6)) 0341 fp_msg (" tiled_gzip_1\n"); 0342 else if (! strncmp (val+1, "GZIP_2", 6)) 0343 fp_msg (" tiled_gzip_2\n"); 0344 else if (! strncmp (val+1, "PLIO_1", 6)) 0345 fp_msg (" tiled_plio\n"); 0346 else if (! strncmp (val+1, "HCOMPRESS_1", 11)) 0347 fp_msg (" tiled_hcompress\n"); 0348 else 0349 fp_msg (" unknown\n"); 0350 0351 } else 0352 fp_msg (" not_tiled\n"); 0353 0354 } else if (hdutype == ASCII_TBL) { 0355 snprintf (msg, SZ_STR," %d ASCII_TBL", hdupos); fp_msg (msg); 0356 snprintf (msg, SZ_STR," SUMS=%lu/%lu\n", (unsigned long) (~((int) hdusum)), datasum); fp_msg (msg); 0357 0358 } else if (hdutype == BINARY_TBL) { 0359 snprintf (msg, SZ_STR," %d BINARY_TBL", hdupos); fp_msg (msg); 0360 snprintf (msg, SZ_STR," SUMS=%lu/%lu\n", (unsigned long) (~((int) hdusum)), datasum); fp_msg (msg); 0361 0362 } else { 0363 snprintf (msg, SZ_STR," %d OTHER", hdupos); fp_msg (msg); 0364 snprintf (msg, SZ_STR," SUMS=%lu/%lu", (unsigned long) (~((int) hdusum)), datasum); fp_msg (msg); 0365 snprintf (msg, SZ_STR," %s\n", val); fp_msg (msg); 0366 } 0367 0368 fits_movrel_hdu (infptr, 1, NULL, &stat); 0369 } 0370 return(0); 0371 } 0372 0373 /*--------------------------------------------------------------------------*/ 0374 int fp_preflight (int argc, char *argv[], int unpack, fpstate *fpptr) 0375 { 0376 char infits[SZ_STR], outfits[SZ_STR]; 0377 int iarg, namelen, nfiles = 0; 0378 0379 if (fpptr->initialized != FP_INIT_MAGIC) { 0380 fp_msg ("Error: internal initialization error\n"); exit (-1); 0381 } 0382 0383 for (iarg=fpptr->firstfile; iarg < argc; iarg++) { 0384 0385 outfits[0] = '\0'; 0386 0387 if (strlen(argv[iarg]) > SZ_STR - 4) { /* allow for .fz or .gz suffix */ 0388 fp_msg ("Error: input file name\n "); fp_msg (argv[iarg]); 0389 fp_msg ("\n is too long\n"); fp_noop (); exit (-1); 0390 } 0391 0392 strncpy (infits, argv[iarg], SZ_STR-1); 0393 if (infits[0] == '-' && infits[1] != '\0') { 0394 /* don't interpret this as intending to read input file from stdin */ 0395 fp_msg ("Error: invalid input file name\n "); fp_msg (argv[iarg]); 0396 fp_msg ("\n"); fp_noop (); exit (-1); 0397 } 0398 0399 if (strchr (infits, '[') || strchr (infits, ']')) { 0400 fp_msg ("Error: section/extension notation not supported: "); 0401 fp_msg (infits); fp_msg ("\n"); fp_noop (); exit (-1); 0402 } 0403 0404 if (unpack) { 0405 /* ********** This section applies to funpack ************ */ 0406 0407 /* check that input file exists */ 0408 if (infits[0] != '-') { /* if not reading from stdin stream */ 0409 if (fp_access (infits) != 0) { /* if not, then check if */ 0410 strcat(infits, ".fz"); /* a .fz version exsits */ 0411 if (fp_access (infits) != 0) { 0412 namelen = strlen(infits); 0413 infits[namelen - 3] = '\0'; /* remove the .fz suffix */ 0414 fp_msg ("Error: can't find or read input file "); fp_msg (infits); 0415 fp_msg ("\n"); fp_noop (); exit (-1); 0416 } 0417 } else { /* make sure a .fz version of the same file doesn't exist */ 0418 namelen = strlen(infits); 0419 strcat(infits, ".fz"); 0420 if (fp_access (infits) == 0) { 0421 infits[namelen] = '\0'; /* remove the .fz suffix */ 0422 fp_msg ("Error: ambiguous input file name. Which file should be unpacked?:\n "); 0423 fp_msg (infits); fp_msg ("\n "); 0424 fp_msg (infits); fp_msg (".fz\n"); 0425 fp_noop (); exit (-1); 0426 } else { 0427 infits[namelen] = '\0'; /* remove the .fz suffix */ 0428 } 0429 } 0430 } 0431 0432 /* if writing to stdout, then we are all done */ 0433 if (fpptr->to_stdout) { 0434 continue; 0435 } 0436 0437 if (fpptr->outfile[0]) { /* user specified output file name */ 0438 nfiles++; 0439 if (nfiles > 1) { 0440 fp_msg ("Error: cannot use same output file name for multiple files:\n "); 0441 fp_msg (fpptr->outfile); 0442 fp_msg ("\n"); fp_noop (); exit (-1); 0443 } 0444 0445 /* check that output file doesn't exist */ 0446 if (fp_access (fpptr->outfile) == 0) { 0447 fp_msg ("Error: output file already exists:\n "); 0448 fp_msg (fpptr->outfile); 0449 fp_msg ("\n "); fp_noop (); exit (-1); 0450 } 0451 continue; 0452 } 0453 0454 /* construct output file name to test */ 0455 if (fpptr->prefix[0]) { 0456 if (strlen(fpptr->prefix) + strlen(infits) > SZ_STR - 1) { 0457 fp_msg ("Error: output file name for\n "); fp_msg (infits); 0458 fp_msg ("\n is too long with the prefix\n"); fp_noop (); exit (-1); 0459 } 0460 strcat(outfits,fpptr->prefix); 0461 } 0462 0463 /* construct output file name */ 0464 if (infits[0] == '-') { 0465 strcpy(outfits, "output.fits"); 0466 } else { 0467 strcpy(outfits, infits); 0468 } 0469 0470 /* remove .gz suffix, if present (output is not gzipped) */ 0471 namelen = strlen(outfits); 0472 if ( !strcmp(".gz", outfits + namelen - 3) ) { 0473 outfits[namelen - 3] = '\0'; 0474 } 0475 0476 /* check for .fz suffix that is sometimes required */ 0477 /* and remove it if present */ 0478 if (infits[0] != '-') { /* if not reading from stdin stream */ 0479 namelen = strlen(outfits); 0480 if ( !strcmp(".fz", outfits + namelen - 3) ) { /* suffix is present */ 0481 outfits[namelen - 3] = '\0'; 0482 } else if (fpptr->delete_suffix) { /* required suffix is missing */ 0483 fp_msg ("Error: input compressed file "); fp_msg (infits); 0484 fp_msg ("\n does not have the default .fz suffix.\n"); 0485 fp_noop (); exit (-1); 0486 } 0487 } 0488 0489 /* if infits != outfits, make sure outfits doesn't already exist */ 0490 if (strcmp(infits, outfits)) { 0491 if (fp_access (outfits) == 0) { 0492 fp_msg ("Error: output file already exists:\n "); fp_msg (outfits); 0493 fp_msg ("\n "); fp_noop (); exit (-1); 0494 } 0495 } 0496 0497 /* if gzipping the output, make sure .gz file doesn't exist */ 0498 if (fpptr->do_gzip_file) { 0499 if (strlen(outfits)+3 > SZ_STR-1) 0500 { 0501 fp_msg ("Error: output file name too long:\n "); fp_msg (outfits); 0502 fp_msg ("\n "); fp_noop (); exit (-1); 0503 } 0504 strcat(outfits, ".gz"); 0505 if (fp_access (outfits) == 0) { 0506 fp_msg ("Error: output file already exists:\n "); fp_msg (outfits); 0507 fp_msg ("\n "); fp_noop (); exit (-1); 0508 } 0509 namelen = strlen(outfits); 0510 outfits[namelen - 3] = '\0'; /* remove the .gz suffix again */ 0511 } 0512 } else { 0513 /* ********** This section applies to fpack ************ */ 0514 0515 /* check that input file exists */ 0516 if (infits[0] != '-') { /* if not reading from stdin stream */ 0517 if (fp_access (infits) != 0) { /* if not, then check if */ 0518 if (strlen(infits)+3 > SZ_STR-1) 0519 { 0520 fp_msg ("Error: input file name too long:\n "); fp_msg (infits); 0521 fp_msg ("\n "); fp_noop (); exit (-1); 0522 } 0523 strcat(infits, ".gz"); /* a gzipped version exsits */ 0524 if (fp_access (infits) != 0) { 0525 namelen = strlen(infits); 0526 infits[namelen - 3] = '\0'; /* remove the .gz suffix */ 0527 fp_msg ("Error: can't find or read input file "); fp_msg (infits); 0528 fp_msg ("\n"); fp_noop (); exit (-1); 0529 } 0530 } 0531 } 0532 0533 /* make sure the file to pack does not already have a .fz suffix */ 0534 namelen = strlen(infits); 0535 if ( !strcmp(".fz", infits + namelen - 3) ) { 0536 fp_msg ("Error: fpack input file already has '.fz' suffix\n" ); fp_msg (infits); 0537 fp_msg ("\n"); fp_noop (); exit (-1); 0538 } 0539 0540 /* if writing to stdout, or just testing the files, then we are all done */ 0541 if (fpptr->to_stdout || fpptr->test_all) { 0542 continue; 0543 } 0544 0545 /* construct output file name */ 0546 if (infits[0] == '-') { 0547 strcpy(outfits, "input.fits"); 0548 } else { 0549 strcpy(outfits, infits); 0550 } 0551 0552 /* remove .gz suffix, if present (output is not gzipped) */ 0553 namelen = strlen(outfits); 0554 if ( !strcmp(".gz", outfits + namelen - 3) ) { 0555 outfits[namelen - 3] = '\0'; 0556 } 0557 0558 /* remove .imh suffix (IRAF format image), and replace with .fits */ 0559 namelen = strlen(outfits); 0560 if ( !strcmp(".imh", outfits + namelen - 4) ) { 0561 outfits[namelen - 4] = '\0'; 0562 strcat(outfits, ".fits"); 0563 } 0564 0565 /* If not clobbering the input file, add .fz suffix to output name */ 0566 if (! fpptr->clobber) 0567 strcat(outfits, ".fz"); 0568 0569 /* if infits != outfits, make sure outfits doesn't already exist */ 0570 if (strcmp(infits, outfits)) { 0571 if (fp_access (outfits) == 0) { 0572 fp_msg ("Error: output file already exists:\n "); fp_msg (outfits); 0573 fp_msg ("\n "); fp_noop (); exit (-1); 0574 } 0575 } 0576 } /* end of fpack section */ 0577 } 0578 0579 fpptr->preflight_checked++; 0580 return(0); 0581 } 0582 0583 /*--------------------------------------------------------------------------*/ 0584 /* must run fp_preflight() before fp_loop() 0585 */ 0586 int fp_loop (int argc, char *argv[], int unpack, char *output_filename, fpstate fpvar) 0587 { 0588 char infits[SZ_STR], outfits[SZ_STR]; 0589 char temp[SZ_STR], answer[30]; 0590 char valchar[]="0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz.#()+,-_@[]/^{}"; 0591 int ichar=0, outlen=0, iarg, islossless, namelen, iraf_infile = 0, status = 0, ifail; 0592 0593 if (fpvar.initialized != FP_INIT_MAGIC) { 0594 fp_msg ("Error: internal initialization error\n"); exit (-1); 0595 } else if (! fpvar.preflight_checked) { 0596 fp_msg ("Error: internal preflight error\n"); exit (-1); 0597 } 0598 0599 if (fpvar.test_all && fpvar.outfile[0]) { 0600 outreport = fopen(fpvar.outfile, "w"); 0601 fprintf(outreport," Filename Extension BITPIX NAXIS1 NAXIS2 Size N_nulls Minval Maxval Mean Sigm Noise1 Noise2 Noise3 Noise5 T_whole T_rowbyrow "); 0602 fprintf(outreport,"[Comp_ratio, Pack_cpu, Unpack_cpu, Lossless readtimes] (repeated for Rice, Hcompress, and GZIP)\n"); 0603 } 0604 0605 0606 tempfilename[0] = '\0'; 0607 tempfilename2[0] = '\0'; 0608 tempfilename3[0] = '\0'; 0609 0610 /* set up signal handler to delete temporary file on abort */ 0611 #ifdef SIGINT 0612 if (signal(SIGINT, SIG_IGN) != SIG_IGN) { 0613 (void) signal(SIGINT, abort_fpack); 0614 } 0615 #endif 0616 0617 #ifdef SIGTERM 0618 if (signal(SIGTERM, SIG_IGN) != SIG_IGN) { 0619 (void) signal(SIGTERM, abort_fpack); 0620 } 0621 #endif 0622 0623 #ifdef SIGHUP 0624 if (signal(SIGHUP, SIG_IGN) != SIG_IGN) { 0625 (void) signal(SIGHUP, abort_fpack); 0626 } 0627 #endif 0628 0629 for (iarg=fpvar.firstfile; iarg < argc; iarg++) { 0630 0631 temp[0] = '\0'; 0632 outfits[0] = '\0'; 0633 islossless = 1; 0634 0635 strncpy (infits, argv[iarg], SZ_STR - 1); 0636 infits[SZ_STR-1]=0; 0637 0638 if (unpack) { 0639 /* ********** This section applies to funpack ************ */ 0640 0641 /* find input file */ 0642 if (infits[0] != '-') { /* if not reading from stdin stream */ 0643 if (fp_access (infits) != 0) { /* if not, then */ 0644 strcat(infits, ".fz"); /* a .fz version must exsit */ 0645 } 0646 } 0647 0648 if (fpvar.to_stdout) { 0649 strcpy(outfits, "-"); 0650 0651 } else if (fpvar.outfile[0]) { /* user specified output file name */ 0652 strcpy(outfits, fpvar.outfile); 0653 0654 } else { 0655 /* construct output file name */ 0656 if (fpvar.prefix[0]) { 0657 strcat(outfits,fpvar.prefix); 0658 } 0659 0660 /* construct output file name */ 0661 if (infits[0] == '-') { 0662 strcpy(outfits, "output.fits"); 0663 } else { 0664 /*strcpy(outfits, infits);*/ 0665 strcpy(outfits, output_filename); 0666 } 0667 0668 /* remove .gz suffix, if present (output is not gzipped) */ 0669 namelen = strlen(outfits); 0670 if ( !strcmp(".gz", outfits + namelen - 3) ) { 0671 outfits[namelen - 3] = '\0'; 0672 } 0673 0674 /* check for .fz suffix that is sometimes required */ 0675 /* and remove it if present */ 0676 namelen = strlen(outfits); 0677 if ( !strcmp(".fz", outfits + namelen - 3) ) { /* suffix is present */ 0678 outfits[namelen - 3] = '\0'; 0679 } 0680 } 0681 0682 } else { 0683 /* ********** This section applies to fpack ************ */ 0684 0685 if (fpvar.to_stdout) { 0686 strcpy(outfits, "-"); 0687 } else if (! fpvar.test_all) { 0688 0689 /* construct output file name */ 0690 if (infits[0] == '-') { 0691 strcpy(outfits, "input.fits"); 0692 } else { 0693 strcpy(outfits, output_filename); 0694 } 0695 0696 /* remove .gz suffix, if present (output is not gzipped) */ 0697 namelen = strlen(outfits); 0698 if ( !strcmp(".gz", outfits + namelen - 3) ) { 0699 outfits[namelen - 3] = '\0'; 0700 } 0701 0702 /* remove .imh suffix (IRAF format image), and replace with .fits */ 0703 namelen = strlen(outfits); 0704 if ( !strcmp(".imh", outfits + namelen - 4) ) { 0705 outfits[namelen - 4] = '\0'; 0706 strcat(outfits, ".fits"); 0707 iraf_infile = 1; /* this is an IRAF format input file */ 0708 /* change the output name to "NAME.fits.fz" */ 0709 } 0710 0711 /* If not clobbering the input file, add .fz suffix to output name */ 0712 if (! fpvar.clobber) 0713 strcat(outfits, ".fz"); 0714 } 0715 } 0716 0717 strncpy(temp, outfits, SZ_STR-1); 0718 temp[SZ_STR-1]=0; 0719 0720 if (infits[0] != '-') { /* if not reading from stdin stream */ 0721 if (!strcmp(infits, outfits) ) { /* are input and output names the same? */ 0722 0723 /* clobber the input file with the output file with the same name */ 0724 if (! fpvar.clobber) { 0725 fp_msg ("\nError: must use -F flag to clobber input file.\n"); 0726 exit (-1); 0727 } 0728 0729 /* create temporary file name in the output directory (same as input directory)*/ 0730 fp_tmpnam("Tmp1", infits, outfits); 0731 0732 strcpy(tempfilename, outfits); /* store temp file name, in case of abort */ 0733 } 0734 } 0735 0736 0737 /* *************** now do the real work ********************* */ 0738 0739 if (fpvar.verbose && ! fpvar.to_stdout) 0740 printf("%s ", infits); 0741 0742 if (fpvar.test_all) { /* compare all the algorithms */ 0743 0744 /* create 2 temporary file names, in the CWD */ 0745 fp_tmpnam("Tmpfile1", "", tempfilename); 0746 fp_tmpnam("Tmpfile2", "", tempfilename2); 0747 0748 fp_test (infits, tempfilename, tempfilename2, fpvar); 0749 0750 remove(tempfilename); 0751 tempfilename[0] = '\0'; /* clear the temp file name */ 0752 remove(tempfilename2); 0753 tempfilename2[0] = '\0'; 0754 continue; 0755 0756 } else if (unpack) { 0757 if (fpvar.to_stdout) { 0758 /* unpack the input file to the stdout stream */ 0759 fp_unpack (infits, outfits, fpvar); 0760 } else { 0761 /* unpack to temporary file, so other tasks can't open it until it is renamed */ 0762 0763 /* create temporary file name, in the output directory */ 0764 fp_tmpnam("Tmp2", outfits, tempfilename2); 0765 0766 /* unpack the input file to the temporary file */ 0767 fp_unpack (infits, tempfilename2, fpvar); 0768 0769 /* rename the temporary file to it's real name */ 0770 ifail = rename(tempfilename2, outfits); 0771 if (ifail) { 0772 fp_msg("Failed to rename temporary file name:\n "); 0773 fp_msg(tempfilename2); 0774 fp_msg(" -> "); 0775 fp_msg(outfits); 0776 fp_msg("\n"); 0777 exit (-1); 0778 } else { 0779 tempfilename2[0] = '\0'; /* clear temporary file name */ 0780 } 0781 } 0782 } else { 0783 fp_pack (infits, outfits, fpvar, &islossless); 0784 } 0785 0786 if (fpvar.to_stdout) { 0787 continue; 0788 } 0789 0790 /* ********** clobber and/or delete files, if needed ************** */ 0791 0792 if (!strcmp(infits, temp) && fpvar.clobber ) { 0793 0794 if (!islossless && ! fpvar.do_not_prompt) { 0795 fp_msg ("\nFile "); 0796 fp_msg (infits); 0797 fp_msg ("\nwas compressed with a LOSSY method. Overwrite the\n"); 0798 fp_msg ("original file with the compressed version? (Y/N) "); 0799 if (fgets(answer, 29, stdin) && answer[0] != 'Y' && answer[0] != 'y') { 0800 fp_msg ("\noriginal file NOT overwritten!\n"); 0801 remove(outfits); 0802 continue; 0803 } 0804 } 0805 0806 if (iraf_infile) { /* special case of deleting an IRAF format header and pixel file */ 0807 if (fits_delete_iraf_file(infits, &status)) { 0808 fp_msg("\nError deleting IRAF .imh and .pix files.\n"); 0809 fp_msg(infits); fp_msg ("\n"); exit (-1); 0810 } 0811 } 0812 0813 #if defined(unix) || defined(__unix__) || defined(__unix) 0814 /* rename clobbers input on Unix platforms */ 0815 if (rename (outfits, temp) != 0) { 0816 fp_msg ("\nError renaming tmp file to "); 0817 fp_msg (temp); fp_msg ("\n"); exit (-1); 0818 } 0819 #else 0820 /* rename DOES NOT clobber existing files on Windows platforms */ 0821 /* so explicitly remove any existing file before renaming the file */ 0822 remove(temp); 0823 if (rename (outfits, temp) != 0) { 0824 fp_msg ("\nError renaming tmp file to "); 0825 fp_msg (temp); fp_msg ("\n"); exit (-1); 0826 } 0827 #endif 0828 0829 tempfilename[0] = '\0'; /* clear temporary file name */ 0830 strcpy(outfits, temp); 0831 0832 } else if (fpvar.clobber || fpvar.delete_input) { /* delete the input file */ 0833 if (!islossless && !fpvar.do_not_prompt) { /* user did not turn off delete prompt */ 0834 fp_msg ("\nFile "); 0835 fp_msg (infits); 0836 fp_msg ("\nwas compressed with a LOSSY method. \n"); 0837 fp_msg ("Delete the original file? (Y/N) "); 0838 if (fgets(answer, 29, stdin) && answer[0] != 'Y' && answer[0] != 'y') { /* user abort */ 0839 fp_msg ("\noriginal file NOT deleted!\n"); 0840 } else { 0841 if (iraf_infile) { /* special case of deleting an IRAF format header and pixel file */ 0842 if (fits_delete_iraf_file(infits, &status)) { 0843 fp_msg("\nError deleting IRAF .imh and .pix files.\n"); 0844 fp_msg(infits); fp_msg ("\n"); exit (-1); 0845 } 0846 } else if (remove(infits) != 0) { /* normal case of deleting input FITS file */ 0847 fp_msg ("\nError deleting input file "); 0848 fp_msg (infits); fp_msg ("\n"); exit (-1); 0849 } 0850 } 0851 } else { /* user said don't prompt, so just delete the input file */ 0852 if (iraf_infile) { /* special case of deleting an IRAF format header and pixel file */ 0853 if (fits_delete_iraf_file(infits, &status)) { 0854 fp_msg("\nError deleting IRAF .imh and .pix files.\n"); 0855 fp_msg(infits); fp_msg ("\n"); exit (-1); 0856 } 0857 } else if (remove(infits) != 0) { /* normal case of deleting input FITS file */ 0858 fp_msg ("\nError deleting input file "); 0859 fp_msg (infits); fp_msg ("\n"); exit (-1); 0860 } 0861 } 0862 } 0863 iraf_infile = 0; 0864 0865 if (fpvar.do_gzip_file) { /* gzip the output file */ 0866 strcpy(temp, "gzip -1 "); 0867 outlen = strlen(outfits); 0868 if (outlen + 8 > SZ_STR-1) 0869 { 0870 fp_msg("\nError: Output file name is too long.\n"); 0871 exit(-1); 0872 } 0873 for (ichar=0; ichar < outlen; ++ichar) 0874 { 0875 if (!strchr(valchar, outfits[ichar])) 0876 { 0877 fp_msg("\n Error: Invalid characters in output file name.\n"); 0878 exit(-1); 0879 } 0880 } 0881 strcat(temp,outfits); 0882 int rc = system(temp); 0883 UNUSED(rc); 0884 strcat(outfits, ".gz"); /* only possibible with funpack */ 0885 } 0886 0887 if (fpvar.verbose && ! fpvar.to_stdout) 0888 printf("-> %s\n", outfits); 0889 0890 } 0891 0892 if (fpvar.test_all && fpvar.outfile[0]) 0893 fclose(outreport); 0894 return(0); 0895 } 0896 0897 /*--------------------------------------------------------------------------*/ 0898 /* fp_pack assumes the output file does not exist (checked by preflight) 0899 */ 0900 int fp_pack (char *infits, char *outfits, fpstate fpvar, int *islossless) 0901 { 0902 fitsfile *infptr, *outfptr; 0903 int stat=0; 0904 0905 fits_open_file (&infptr, infits, READONLY, &stat); 0906 if (stat) 0907 { 0908 fits_report_error (stderr, stat); 0909 return -1; 0910 } 0911 0912 fits_create_file (&outfptr, outfits, &stat); 0913 if (stat) 0914 { 0915 fp_abort_output(infptr, outfptr, stat); 0916 return -1; 0917 } 0918 0919 while (! stat) { 0920 0921 /* LOOP OVER EACH HDU */ 0922 0923 fits_set_lossy_int (outfptr, fpvar.int_to_float, &stat); 0924 fits_set_compression_type (outfptr, fpvar.comptype, &stat); 0925 fits_set_tile_dim (outfptr, 6, fpvar.ntile, &stat); 0926 0927 if (fpvar.no_dither) 0928 fits_set_quantize_method(outfptr, -1, &stat); 0929 else 0930 fits_set_quantize_method(outfptr, fpvar.dither_method, &stat); 0931 0932 fits_set_quantize_level (outfptr, fpvar.quantize_level, &stat); 0933 fits_set_dither_offset(outfptr, fpvar.dither_offset, &stat); 0934 fits_set_hcomp_scale (outfptr, fpvar.scale, &stat); 0935 fits_set_hcomp_smooth (outfptr, fpvar.smooth, &stat); 0936 0937 fp_pack_hdu (infptr, outfptr, fpvar, islossless, &stat); 0938 0939 if (fpvar.do_checksums) { 0940 fits_write_chksum (outfptr, &stat); 0941 } 0942 0943 fits_movrel_hdu (infptr, 1, NULL, &stat); 0944 } 0945 0946 if (stat == END_OF_FILE) stat = 0; 0947 0948 /* set checksum for case of newly created primary HDU */ 0949 0950 if (fpvar.do_checksums) { 0951 fits_movabs_hdu (outfptr, 1, NULL, &stat); 0952 fits_write_chksum (outfptr, &stat); 0953 } 0954 0955 if (stat) { 0956 fp_abort_output(infptr, outfptr, stat); 0957 } 0958 0959 fits_close_file (outfptr, &stat); 0960 fits_close_file (infptr, &stat); 0961 0962 return(0); 0963 } 0964 0965 /*--------------------------------------------------------------------------*/ 0966 0967 int fp_pack_fits_to_fits (fitsfile *infptr, fitsfile **outfits, fpstate fpvar, int *islossless) 0968 { 0969 fitsfile *outfptr; 0970 int stat=0; 0971 size_t outbufferSize = 2880; 0972 void *outbuffer = (char *)malloc(outbufferSize); 0973 0974 fits_create_memfile(&outfptr, &outbuffer, &outbufferSize, 2880, realloc, &stat); 0975 if (stat) 0976 { 0977 fp_abort_output(infptr, outfptr, stat); 0978 return -1; 0979 } 0980 0981 while (! stat) { 0982 0983 /* LOOP OVER EACH HDU */ 0984 0985 fits_set_lossy_int (outfptr, fpvar.int_to_float, &stat); 0986 fits_set_compression_type (outfptr, fpvar.comptype, &stat); 0987 fits_set_tile_dim (outfptr, 6, fpvar.ntile, &stat); 0988 0989 if (fpvar.no_dither) 0990 fits_set_quantize_method(outfptr, -1, &stat); 0991 else 0992 fits_set_quantize_method(outfptr, fpvar.dither_method, &stat); 0993 0994 fits_set_quantize_level (outfptr, fpvar.quantize_level, &stat); 0995 fits_set_dither_offset(outfptr, fpvar.dither_offset, &stat); 0996 fits_set_hcomp_scale (outfptr, fpvar.scale, &stat); 0997 fits_set_hcomp_smooth (outfptr, fpvar.smooth, &stat); 0998 0999 fp_pack_hdu (infptr, outfptr, fpvar, islossless, &stat); 1000 1001 if (fpvar.do_checksums) { 1002 fits_write_chksum (outfptr, &stat); 1003 } 1004 1005 fits_movrel_hdu (infptr, 1, NULL, &stat); 1006 } 1007 1008 if (stat == END_OF_FILE) stat = 0; 1009 1010 /* set checksum for case of newly created primary HDU */ 1011 1012 if (fpvar.do_checksums) { 1013 fits_movabs_hdu (outfptr, 1, NULL, &stat); 1014 fits_write_chksum (outfptr, &stat); 1015 } 1016 1017 if (stat) { 1018 fp_abort_output(infptr, outfptr, stat); 1019 } 1020 1021 fits_close_file (infptr, &stat); 1022 *outfits = outfptr; 1023 1024 return(0); 1025 } 1026 1027 /*--------------------------------------------------------------------------*/ 1028 /* fp_unpack assumes the output file does not exist 1029 */ 1030 int fp_unpack (char *infits, char *outfits, fpstate fpvar) 1031 { 1032 fitsfile *infptr, *outfptr; 1033 int stat=0, hdutype, extnum, single = 0; 1034 char *loc, *hduloc = 0, hduname[SZ_STR] = { 0 }; 1035 1036 fits_open_file (&infptr, infits, READONLY, &stat); 1037 fits_create_file (&outfptr, outfits, &stat); 1038 1039 if (stat) { 1040 fp_abort_output(infptr, outfptr, stat); 1041 } 1042 1043 if (fpvar.extname[0]) { /* unpack a list of HDUs? */ 1044 1045 /* move to the first HDU in the list */ 1046 hduloc = fpvar.extname; 1047 loc = strchr(hduloc, ','); /* look for 'comma' delimiter between names */ 1048 1049 if (loc) 1050 *loc = '\0'; /* terminate the first name in the string */ 1051 1052 strcpy(hduname, hduloc); /* copy the first name into temporary string */ 1053 1054 if (loc) 1055 hduloc = loc + 1; /* advance to the beginning of the next name, if any */ 1056 else { 1057 hduloc += strlen(hduname); /* end of the list */ 1058 single = 1; /* only 1 HDU is being unpacked */ 1059 } 1060 1061 if (isdigit( (int) hduname[0]) ) { 1062 extnum = strtol(hduname, &loc, 10); /* read the string as an integer */ 1063 1064 /* check for junk following the integer */ 1065 if (*loc == '\0' ) /* no junk, so move to this HDU number (+1) */ 1066 { 1067 fits_movabs_hdu(infptr, extnum + 1, &hdutype, &stat); /* move to HDU number */ 1068 if (hdutype != IMAGE_HDU) 1069 stat = NOT_IMAGE; 1070 1071 } else { /* the string is not an integer, so must be the column name */ 1072 hdutype = IMAGE_HDU; 1073 fits_movnam_hdu(infptr, hdutype, hduname, 0, &stat); 1074 } 1075 } 1076 else 1077 { 1078 /* move to the named image extension */ 1079 hdutype = IMAGE_HDU; 1080 fits_movnam_hdu(infptr, hdutype, hduname, 0, &stat); 1081 } 1082 } 1083 1084 if (stat) { 1085 fp_msg ("Unable to find and move to extension '"); 1086 fp_msg(hduname); 1087 fp_msg("'\n"); 1088 fp_abort_output(infptr, outfptr, stat); 1089 } 1090 1091 while (! stat) { 1092 1093 if (single) 1094 stat = -1; /* special status flag to force output primary array */ 1095 1096 fp_unpack_hdu (infptr, outfptr, fpvar, &stat); 1097 1098 if (fpvar.do_checksums) { 1099 fits_write_chksum (outfptr, &stat); 1100 } 1101 1102 /* move to the next HDU */ 1103 if (fpvar.extname[0]) { /* unpack a list of HDUs? */ 1104 1105 if (!(*hduloc)) { 1106 stat = END_OF_FILE; /* we reached the end of the list */ 1107 } else { 1108 /* parse the next HDU name and move to it */ 1109 loc = strchr(hduloc, ','); 1110 1111 if (loc) /* look for 'comma' delimiter between names */ 1112 *loc = '\0'; /* terminate the first name in the string */ 1113 1114 strcpy(hduname, hduloc); /* copy the next name into temporary string */ 1115 1116 if (loc) 1117 hduloc = loc + 1; /* advance to the beginning of the next name, if any */ 1118 else 1119 *hduloc = '\0'; /* end of the list */ 1120 1121 if (isdigit( (int) hduname[0]) ) { 1122 extnum = strtol(hduname, &loc, 10); /* read the string as an integer */ 1123 1124 /* check for junk following the integer */ 1125 if (*loc == '\0' ) /* no junk, so move to this HDU number (+1) */ 1126 { 1127 fits_movabs_hdu(infptr, extnum + 1, &hdutype, &stat); /* move to HDU number */ 1128 if (hdutype != IMAGE_HDU) 1129 stat = NOT_IMAGE; 1130 1131 } else { /* the string is not an integer, so must be the column name */ 1132 hdutype = IMAGE_HDU; 1133 fits_movnam_hdu(infptr, hdutype, hduname, 0, &stat); 1134 } 1135 1136 } else { 1137 /* move to the named image extension */ 1138 hdutype = IMAGE_HDU; 1139 fits_movnam_hdu(infptr, hdutype, hduname, 0, &stat); 1140 } 1141 1142 if (stat) { 1143 fp_msg ("Unable to find and move to extension '"); 1144 fp_msg(hduname); 1145 fp_msg("'\n"); 1146 } 1147 } 1148 } else { 1149 /* increment to the next HDU */ 1150 fits_movrel_hdu (infptr, 1, NULL, &stat); 1151 } 1152 } 1153 1154 if (stat == END_OF_FILE) stat = 0; 1155 1156 /* set checksum for case of newly created primary HDU 1157 */ 1158 if (fpvar.do_checksums) { 1159 fits_movabs_hdu (outfptr, 1, NULL, &stat); 1160 fits_write_chksum (outfptr, &stat); 1161 } 1162 1163 1164 if (stat) { 1165 fp_abort_output(infptr, outfptr, stat); 1166 } 1167 1168 fits_close_file (outfptr, &stat); 1169 fits_close_file (infptr, &stat); 1170 1171 return(0); 1172 } 1173 1174 /*--------------------------------------------------------------------------*/ 1175 /* fp_unpack assumes the output file does not exist 1176 */ 1177 int fp_unpack_file_to_fits (char *infits, fitsfile **outfits, fpstate fpvar) 1178 { 1179 fitsfile *infptr, *outfptr; 1180 int stat=0, hdutype, extnum, single = 0; 1181 char *loc, *hduloc = 0, hduname[SZ_STR] = { 0 }; 1182 size_t outbufferSize = 2880; 1183 void *outbuffer = (char *)malloc(outbufferSize); 1184 1185 fits_open_file (&infptr, infits, READONLY, &stat); 1186 fits_create_memfile(&outfptr, &outbuffer, &outbufferSize, 2880, realloc, &stat); 1187 1188 if (stat) { 1189 fp_abort_output(infptr, outfptr, stat); 1190 } 1191 1192 if (fpvar.extname[0]) { /* unpack a list of HDUs? */ 1193 1194 /* move to the first HDU in the list */ 1195 hduloc = fpvar.extname; 1196 loc = strchr(hduloc, ','); /* look for 'comma' delimiter between names */ 1197 1198 if (loc) 1199 *loc = '\0'; /* terminate the first name in the string */ 1200 1201 strcpy(hduname, hduloc); /* copy the first name into temporary string */ 1202 1203 if (loc) 1204 hduloc = loc + 1; /* advance to the beginning of the next name, if any */ 1205 else { 1206 hduloc += strlen(hduname); /* end of the list */ 1207 single = 1; /* only 1 HDU is being unpacked */ 1208 } 1209 1210 if (isdigit( (int) hduname[0]) ) { 1211 extnum = strtol(hduname, &loc, 10); /* read the string as an integer */ 1212 1213 /* check for junk following the integer */ 1214 if (*loc == '\0' ) /* no junk, so move to this HDU number (+1) */ 1215 { 1216 fits_movabs_hdu(infptr, extnum + 1, &hdutype, &stat); /* move to HDU number */ 1217 if (hdutype != IMAGE_HDU) 1218 stat = NOT_IMAGE; 1219 1220 } else { /* the string is not an integer, so must be the column name */ 1221 hdutype = IMAGE_HDU; 1222 fits_movnam_hdu(infptr, hdutype, hduname, 0, &stat); 1223 } 1224 } 1225 else 1226 { 1227 /* move to the named image extension */ 1228 hdutype = IMAGE_HDU; 1229 fits_movnam_hdu(infptr, hdutype, hduname, 0, &stat); 1230 } 1231 } 1232 1233 if (stat) { 1234 fp_msg ("Unable to find and move to extension '"); 1235 fp_msg(hduname); 1236 fp_msg("'\n"); 1237 fp_abort_output(infptr, outfptr, stat); 1238 } 1239 1240 while (! stat) { 1241 1242 if (single) 1243 stat = -1; /* special status flag to force output primary array */ 1244 1245 fp_unpack_hdu (infptr, outfptr, fpvar, &stat); 1246 1247 if (fpvar.do_checksums) { 1248 fits_write_chksum (outfptr, &stat); 1249 } 1250 1251 /* move to the next HDU */ 1252 if (fpvar.extname[0]) { /* unpack a list of HDUs? */ 1253 1254 if (!(*hduloc)) { 1255 stat = END_OF_FILE; /* we reached the end of the list */ 1256 } else { 1257 /* parse the next HDU name and move to it */ 1258 loc = strchr(hduloc, ','); 1259 1260 if (loc) /* look for 'comma' delimiter between names */ 1261 *loc = '\0'; /* terminate the first name in the string */ 1262 1263 strcpy(hduname, hduloc); /* copy the next name into temporary string */ 1264 1265 if (loc) 1266 hduloc = loc + 1; /* advance to the beginning of the next name, if any */ 1267 else 1268 *hduloc = '\0'; /* end of the list */ 1269 1270 if (isdigit( (int) hduname[0]) ) { 1271 extnum = strtol(hduname, &loc, 10); /* read the string as an integer */ 1272 1273 /* check for junk following the integer */ 1274 if (*loc == '\0' ) /* no junk, so move to this HDU number (+1) */ 1275 { 1276 fits_movabs_hdu(infptr, extnum + 1, &hdutype, &stat); /* move to HDU number */ 1277 if (hdutype != IMAGE_HDU) 1278 stat = NOT_IMAGE; 1279 1280 } else { /* the string is not an integer, so must be the column name */ 1281 hdutype = IMAGE_HDU; 1282 fits_movnam_hdu(infptr, hdutype, hduname, 0, &stat); 1283 } 1284 1285 } else { 1286 /* move to the named image extension */ 1287 hdutype = IMAGE_HDU; 1288 fits_movnam_hdu(infptr, hdutype, hduname, 0, &stat); 1289 } 1290 1291 if (stat) { 1292 fp_msg ("Unable to find and move to extension '"); 1293 fp_msg(hduname); 1294 fp_msg("'\n"); 1295 } 1296 } 1297 } else { 1298 /* increment to the next HDU */ 1299 fits_movrel_hdu (infptr, 1, NULL, &stat); 1300 } 1301 } 1302 1303 if (stat == END_OF_FILE) stat = 0; 1304 1305 /* set checksum for case of newly created primary HDU 1306 */ 1307 if (fpvar.do_checksums) { 1308 fits_movabs_hdu (outfptr, 1, NULL, &stat); 1309 fits_write_chksum (outfptr, &stat); 1310 } 1311 1312 1313 if (stat) { 1314 fp_abort_output(infptr, outfptr, stat); 1315 } 1316 1317 fits_close_file (infptr, &stat); 1318 *outfits = outfptr; 1319 1320 return(0); 1321 } 1322 1323 /*--------------------------------------------------------------------------*/ 1324 /* fp_unpack_memroy is like fp_unpack but unpacks in RAM. 1325 */ 1326 int fp_unpack_data_to_data (const char *inputBuffer, size_t inputBufferSize, unsigned char **outputBuffer, 1327 size_t *outputBufferSize, fpstate fpvar) 1328 { 1329 fitsfile *infptr, *outfptr; 1330 int stat=0, hdutype, extnum, single = 0; 1331 char *loc, *hduloc = 0, hduname[SZ_STR] = { 0 }; 1332 1333 fits_open_memfile(&infptr, "", READONLY, &inputBuffer, &inputBufferSize, 2880, NULL, &stat); 1334 fits_create_memfile(&outfptr, outputBuffer, outputBufferSize, 100000, realloc, &stat); 1335 1336 if (stat) { 1337 fp_abort_output(infptr, outfptr, stat); 1338 } 1339 1340 if (fpvar.extname[0]) { /* unpack a list of HDUs? */ 1341 1342 /* move to the first HDU in the list */ 1343 hduloc = fpvar.extname; 1344 loc = strchr(hduloc, ','); /* look for 'comma' delimiter between names */ 1345 1346 if (loc) 1347 *loc = '\0'; /* terminate the first name in the string */ 1348 1349 strcpy(hduname, hduloc); /* copy the first name into temporary string */ 1350 1351 if (loc) 1352 hduloc = loc + 1; /* advance to the beginning of the next name, if any */ 1353 else { 1354 hduloc += strlen(hduname); /* end of the list */ 1355 single = 1; /* only 1 HDU is being unpacked */ 1356 } 1357 1358 if (isdigit( (int) hduname[0]) ) { 1359 extnum = strtol(hduname, &loc, 10); /* read the string as an integer */ 1360 1361 /* check for junk following the integer */ 1362 if (*loc == '\0' ) /* no junk, so move to this HDU number (+1) */ 1363 { 1364 fits_movabs_hdu(infptr, extnum + 1, &hdutype, &stat); /* move to HDU number */ 1365 if (hdutype != IMAGE_HDU) 1366 stat = NOT_IMAGE; 1367 1368 } else { /* the string is not an integer, so must be the column name */ 1369 hdutype = IMAGE_HDU; 1370 fits_movnam_hdu(infptr, hdutype, hduname, 0, &stat); 1371 } 1372 } 1373 else 1374 { 1375 /* move to the named image extension */ 1376 hdutype = IMAGE_HDU; 1377 fits_movnam_hdu(infptr, hdutype, hduname, 0, &stat); 1378 } 1379 } 1380 1381 if (stat) { 1382 fp_msg ("Unable to find and move to extension '"); 1383 fp_msg(hduname); 1384 fp_msg("'\n"); 1385 fp_abort_output(infptr, outfptr, stat); 1386 } 1387 1388 while (! stat) { 1389 1390 if (single) 1391 stat = -1; /* special status flag to force output primary array */ 1392 1393 fp_unpack_hdu (infptr, outfptr, fpvar, &stat); 1394 1395 if (fpvar.do_checksums) { 1396 fits_write_chksum (outfptr, &stat); 1397 } 1398 1399 /* move to the next HDU */ 1400 if (fpvar.extname[0]) { /* unpack a list of HDUs? */ 1401 1402 if (!(*hduloc)) { 1403 stat = END_OF_FILE; /* we reached the end of the list */ 1404 } else { 1405 /* parse the next HDU name and move to it */ 1406 loc = strchr(hduloc, ','); 1407 1408 if (loc) /* look for 'comma' delimiter between names */ 1409 *loc = '\0'; /* terminate the first name in the string */ 1410 1411 strcpy(hduname, hduloc); /* copy the next name into temporary string */ 1412 1413 if (loc) 1414 hduloc = loc + 1; /* advance to the beginning of the next name, if any */ 1415 else 1416 *hduloc = '\0'; /* end of the list */ 1417 1418 if (isdigit( (int) hduname[0]) ) { 1419 extnum = strtol(hduname, &loc, 10); /* read the string as an integer */ 1420 1421 /* check for junk following the integer */ 1422 if (*loc == '\0' ) /* no junk, so move to this HDU number (+1) */ 1423 { 1424 fits_movabs_hdu(infptr, extnum + 1, &hdutype, &stat); /* move to HDU number */ 1425 if (hdutype != IMAGE_HDU) 1426 stat = NOT_IMAGE; 1427 1428 } else { /* the string is not an integer, so must be the column name */ 1429 hdutype = IMAGE_HDU; 1430 fits_movnam_hdu(infptr, hdutype, hduname, 0, &stat); 1431 } 1432 1433 } else { 1434 /* move to the named image extension */ 1435 hdutype = IMAGE_HDU; 1436 fits_movnam_hdu(infptr, hdutype, hduname, 0, &stat); 1437 } 1438 1439 if (stat) { 1440 fp_msg ("Unable to find and move to extension '"); 1441 fp_msg(hduname); 1442 fp_msg("'\n"); 1443 } 1444 } 1445 } else { 1446 /* increment to the next HDU */ 1447 fits_movrel_hdu (infptr, 1, NULL, &stat); 1448 } 1449 } 1450 1451 if (stat == END_OF_FILE) stat = 0; 1452 1453 /* set checksum for case of newly created primary HDU 1454 */ 1455 if (fpvar.do_checksums) { 1456 fits_movabs_hdu (outfptr, 1, NULL, &stat); 1457 fits_write_chksum (outfptr, &stat); 1458 } 1459 1460 if (stat) { 1461 fp_abort_output(infptr, outfptr, stat); 1462 } 1463 1464 fits_close_file (outfptr, &stat); 1465 fits_close_file (infptr, &stat); 1466 return(0); 1467 } 1468 1469 /*--------------------------------------------------------------------------*/ 1470 /* fp_test assumes the output files do not exist 1471 */ 1472 int fp_test (char *infits, char *outfits, char *outfits2, fpstate fpvar) 1473 { 1474 fitsfile *inputfptr, *infptr, *outfptr, *outfptr2, *tempfile; 1475 1476 long naxes[9] = {1, 1, 1, 1, 1, 1, 1, 1, 1}; 1477 int stat=0, totpix=0, naxis=0, ii, hdutype, bitpix = 0, extnum = 0, len; 1478 int tstatus = 0, hdunum, rescale_flag, bpix = 8, ncols; 1479 char dtype[8], dimen[100]; 1480 double bscale, rescale, noisemin; 1481 long headstart, datastart, dataend; 1482 float origdata = 0., whole_cpu, whole_elapse, row_elapse, row_cpu, xbits; 1483 1484 LONGLONG nrows; 1485 /* structure to hold image statistics (defined in fpack.h) */ 1486 imgstats imagestats = { 0 }; 1487 1488 fits_open_file (&inputfptr, infits, READONLY, &stat); 1489 fits_create_file (&outfptr, outfits, &stat); 1490 fits_create_file (&outfptr2, outfits2, &stat); 1491 1492 if (stat) { fits_report_error (stderr, stat); exit (stat); } 1493 1494 while (! stat) { 1495 1496 /* LOOP OVER EACH HDU */ 1497 rescale_flag = 0; 1498 fits_get_hdu_type (inputfptr, &hdutype, &stat); 1499 1500 if (hdutype == IMAGE_HDU) { 1501 fits_get_img_param (inputfptr, 9, &bitpix, &naxis, naxes, &stat); 1502 for (totpix=1, ii=0; ii < 9; ii++) totpix *= naxes[ii]; 1503 } 1504 1505 if (!fits_is_compressed_image (inputfptr, &stat) && hdutype == IMAGE_HDU && 1506 naxis != 0 && totpix != 0 && fpvar.do_images) { 1507 1508 /* rescale a scaled integer image to reduce noise? */ 1509 if (fpvar.rescale_noise != 0. && bitpix > 0 && bitpix < LONGLONG_IMG) { 1510 1511 tstatus = 0; 1512 fits_read_key(inputfptr, TDOUBLE, "BSCALE", &bscale, 0, &tstatus); 1513 1514 if (tstatus == 0 && bscale != 1.0) { /* image must be scaled */ 1515 1516 if (bitpix == LONG_IMG) 1517 fp_i4stat(inputfptr, naxis, naxes, &imagestats, &stat); 1518 else 1519 fp_i2stat(inputfptr, naxis, naxes, &imagestats, &stat); 1520 1521 /* use the minimum of the MAD 2nd, 3rd, and 5th order noise estimates */ 1522 noisemin = imagestats.noise3; 1523 if (imagestats.noise2 != 0. && imagestats.noise2 < noisemin) noisemin = imagestats.noise2; 1524 if (imagestats.noise5 != 0. && imagestats.noise5 < noisemin) noisemin = imagestats.noise5; 1525 1526 rescale = noisemin / fpvar.rescale_noise; 1527 if (rescale > 1.0) { 1528 1529 /* all the criteria are met, so create a temporary file that */ 1530 /* contains a rescaled version of the image, in CWD */ 1531 1532 /* create temporary file name */ 1533 fp_tmpnam("Tmpfile3", "", tempfilename3); 1534 1535 fits_create_file(&tempfile, tempfilename3, &stat); 1536 1537 fits_get_hdu_num(inputfptr, &hdunum); 1538 if (hdunum != 1) { 1539 1540 /* the input hdu is an image extension, so create dummy primary */ 1541 fits_create_img(tempfile, 8, 0, naxes, &stat); 1542 } 1543 1544 fits_copy_header(inputfptr, tempfile, &stat); /* copy the header */ 1545 1546 /* rescale the data, so that it will compress more efficiently */ 1547 if (bitpix == LONG_IMG) 1548 fp_i4rescale(inputfptr, naxis, naxes, rescale, tempfile, &stat); 1549 else 1550 fp_i2rescale(inputfptr, naxis, naxes, rescale, tempfile, &stat); 1551 1552 /* scale the BSCALE keyword by the inverse factor */ 1553 1554 bscale = bscale * rescale; 1555 fits_update_key(tempfile, TDOUBLE, "BSCALE", &bscale, 0, &stat); 1556 1557 /* rescan the header, to reset the actual scaling parameters */ 1558 fits_set_hdustruc(tempfile, &stat); 1559 1560 infptr = tempfile; 1561 rescale_flag = 1; 1562 } 1563 } 1564 } 1565 1566 if (!rescale_flag) /* just compress the input file, without rescaling */ 1567 infptr = inputfptr; 1568 1569 /* compute basic statistics about the input image */ 1570 if (bitpix == BYTE_IMG) { 1571 bpix = 8; 1572 strncpy(dtype, "8 ", sizeof(dtype)/sizeof(dtype[0])); 1573 fp_i2stat(infptr, naxis, naxes, &imagestats, &stat); 1574 } else if (bitpix == SHORT_IMG) { 1575 bpix = 16; 1576 strncpy(dtype, "16 ", sizeof(dtype)/sizeof(dtype[0])); 1577 fp_i2stat(infptr, naxis, naxes, &imagestats, &stat); 1578 } else if (bitpix == LONG_IMG) { 1579 bpix = 32; 1580 strncpy(dtype, "32 ", sizeof(dtype)/sizeof(dtype[0])); 1581 fp_i4stat(infptr, naxis, naxes, &imagestats, &stat); 1582 } else if (bitpix == LONGLONG_IMG) { 1583 bpix = 64; 1584 strncpy(dtype, "64 ", sizeof(dtype)/sizeof(dtype[0])); 1585 } else if (bitpix == FLOAT_IMG) { 1586 bpix = 32; 1587 strncpy(dtype, "-32", sizeof(dtype)/sizeof(dtype[0])); 1588 fp_r4stat(infptr, naxis, naxes, &imagestats, &stat); 1589 } else if (bitpix == DOUBLE_IMG) { 1590 bpix = 64; 1591 strncpy(dtype, "-64", sizeof(dtype)/sizeof(dtype[0])); 1592 fp_r4stat(infptr, naxis, naxes, &imagestats, &stat); 1593 } 1594 else dtype[0] = '\0'; 1595 1596 /* use the minimum of the MAD 2nd, 3rd, and 5th order noise estimates */ 1597 noisemin = imagestats.noise3; 1598 if (imagestats.noise2 != 0. && imagestats.noise2 < noisemin) noisemin = imagestats.noise2; 1599 if (imagestats.noise5 != 0. && imagestats.noise5 < noisemin) noisemin = imagestats.noise5; 1600 1601 xbits = (float) (log10(noisemin)/.301 + 1.792); 1602 1603 printf("\n File: %s\n", infits); 1604 printf(" Ext BITPIX Dimens. Nulls Min Max Mean Sigma Noise2 Noise3 Noise5 Nbits MaxR\n"); 1605 1606 printf(" %3d %s", extnum, dtype); 1607 snprintf(dimen,100," (%ld", naxes[0]); 1608 len =strlen(dimen); 1609 for (ii = 1; ii < naxis; ii++) { 1610 if (len < 99) 1611 snprintf(dimen+len,100-len,",%ld", naxes[ii]); 1612 len =strlen(dimen); 1613 } 1614 if (strlen(dimen)<99) 1615 strcat(dimen, ")"); 1616 printf("%-12s",dimen); 1617 1618 fits_get_hduaddr(inputfptr, &headstart, &datastart, &dataend, &stat); 1619 origdata = (float) ((dataend - datastart)/1000000.); 1620 1621 /* get elapsed and cpu times need to read the uncompressed image */ 1622 fits_read_image_speed (infptr, &whole_elapse, &whole_cpu, 1623 &row_elapse, &row_cpu, &stat); 1624 1625 printf(" %5d %6.0f %6.0f %8.1f %#8.2g %#7.3g %#7.3g %#7.3g %#5.1f %#6.2f\n", 1626 imagestats.n_nulls, imagestats.minval, imagestats.maxval, 1627 imagestats.mean, imagestats.sigma, 1628 imagestats.noise2, imagestats.noise3, imagestats.noise5, xbits, bpix/xbits); 1629 1630 printf("\n Type Ratio Size (MB) Pk (Sec) UnPk Exact ElpN CPUN Elp1 CPU1\n"); 1631 1632 printf(" Native %5.3f %5.3f %5.3f %5.3f\n", 1633 whole_elapse, whole_cpu, row_elapse, row_cpu); 1634 1635 if (fpvar.outfile[0]) { 1636 fprintf(outreport, 1637 " %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", 1638 infits, extnum, bitpix, naxes[0], naxes[1], origdata, imagestats.n_nulls, imagestats.minval, 1639 imagestats.maxval, imagestats.mean, imagestats.sigma, 1640 imagestats.noise1, imagestats.noise2, imagestats.noise3, imagestats.noise5, whole_elapse, whole_cpu, row_elapse, row_cpu); 1641 } 1642 1643 fits_set_lossy_int (outfptr, fpvar.int_to_float, &stat); 1644 if ( (bitpix > 0) && (fpvar.int_to_float != 0) ) { 1645 1646 if ( (noisemin < (fpvar.n3ratio * fpvar.quantize_level) ) || 1647 (noisemin < fpvar.n3min)) { 1648 1649 /* image contains too little noise to quantize effectively */ 1650 fits_set_lossy_int (outfptr, 0, &stat); 1651 fits_get_hdu_num(infptr, &hdunum); 1652 1653 printf(" HDU %d does not meet noise criteria to be quantized, so losslessly compressed.\n", hdunum); 1654 } 1655 } 1656 1657 /* test compression ratio and speed for each algorithm */ 1658 1659 if (fpvar.quantize_level != 0) { 1660 1661 fits_set_compression_type (outfptr, RICE_1, &stat); 1662 fits_set_tile_dim (outfptr, 6, fpvar.ntile, &stat); 1663 if (fpvar.no_dither) 1664 fits_set_quantize_method(outfptr, -1, &stat); 1665 else 1666 fits_set_quantize_method(outfptr, fpvar.dither_method, &stat); 1667 1668 fits_set_quantize_level (outfptr, fpvar.quantize_level, &stat); 1669 fits_set_dither_offset(outfptr, fpvar.dither_offset, &stat); 1670 fits_set_hcomp_scale (outfptr, fpvar.scale, &stat); 1671 fits_set_hcomp_smooth (outfptr, fpvar.smooth, &stat); 1672 1673 fp_test_hdu(infptr, outfptr, outfptr2, fpvar, &stat); 1674 } 1675 1676 if (fpvar.quantize_level != 0) { 1677 \ 1678 fits_set_compression_type (outfptr, HCOMPRESS_1, &stat); 1679 fits_set_tile_dim (outfptr, 6, fpvar.ntile, &stat); 1680 1681 if (fpvar.no_dither) 1682 fits_set_quantize_method(outfptr, -1, &stat); 1683 else 1684 fits_set_quantize_method(outfptr, fpvar.dither_method, &stat); 1685 1686 fits_set_quantize_level (outfptr, fpvar.quantize_level, &stat); 1687 fits_set_dither_offset(outfptr, fpvar.dither_offset, &stat); 1688 fits_set_hcomp_scale (outfptr, fpvar.scale, &stat); 1689 fits_set_hcomp_smooth (outfptr, fpvar.smooth, &stat); 1690 1691 fp_test_hdu(infptr, outfptr, outfptr2, fpvar, &stat); 1692 } 1693 1694 if (fpvar.comptype == GZIP_2) { 1695 fits_set_compression_type (outfptr, GZIP_2, &stat); 1696 } else { 1697 fits_set_compression_type (outfptr, GZIP_1, &stat); 1698 } 1699 1700 fits_set_tile_dim (outfptr, 6, fpvar.ntile, &stat); 1701 1702 if (fpvar.no_dither) 1703 fits_set_quantize_method(outfptr, -1, &stat); 1704 else 1705 fits_set_quantize_method(outfptr, fpvar.dither_method, &stat); 1706 1707 fits_set_quantize_level (outfptr, fpvar.quantize_level, &stat); 1708 fits_set_dither_offset(outfptr, fpvar.dither_offset, &stat); 1709 fits_set_hcomp_scale (outfptr, fpvar.scale, &stat); 1710 fits_set_hcomp_smooth (outfptr, fpvar.smooth, &stat); 1711 1712 fp_test_hdu(infptr, outfptr, outfptr2, fpvar, &stat); 1713 1714 /* 1715 fits_set_compression_type (outfptr, BZIP2_1, &stat); 1716 fits_set_tile_dim (outfptr, 6, fpvar.ntile, &stat); 1717 fp_test_hdu(infptr, outfptr, outfptr2, fpvar, &stat); 1718 */ 1719 /* 1720 fits_set_compression_type (outfptr, PLIO_1, &stat); 1721 fits_set_tile_dim (outfptr, 6, fpvar.ntile, &stat); 1722 fp_test_hdu(infptr, outfptr, outfptr2, fpvar, &stat); 1723 */ 1724 /* 1725 if (bitpix == SHORT_IMG || bitpix == LONG_IMG) { 1726 fits_set_compression_type (outfptr, NOCOMPRESS, &stat); 1727 fits_set_tile_dim (outfptr, 6, fpvar.ntile, &stat); 1728 fp_test_hdu(infptr, outfptr, outfptr2, fpvar, &stat); 1729 } 1730 */ 1731 if (fpvar.outfile[0]) 1732 fprintf(outreport,"\n"); 1733 1734 /* delete the temporary file */ 1735 if (rescale_flag) { 1736 fits_delete_file (infptr, &stat); 1737 tempfilename3[0] = '\0'; /* clear the temp filename */ 1738 } 1739 } else if ( (hdutype == BINARY_TBL) && fpvar.do_tables) { 1740 1741 fits_get_num_rowsll(inputfptr, &nrows, &stat); 1742 fits_get_num_cols(inputfptr, &ncols, &stat); 1743 #if defined(_MSC_VER) 1744 /* Microsoft Visual C++ 6.0 uses '%I64d' syntax for 8-byte integers */ 1745 printf("\n File: %s, HDU %d, %d cols X %I64d rows\n", infits, extnum, ncols, nrows); 1746 #elif (USE_LL_SUFFIX == 1) 1747 printf("\n File: %s, HDU %d, %d cols X %lld rows\n", infits, extnum, ncols, nrows); 1748 #else 1749 printf("\n File: %s, HDU %d, %d cols X %ld rows\n", infits, extnum, ncols, nrows); 1750 #endif 1751 fp_test_table(inputfptr, outfptr, outfptr2, fpvar, &stat); 1752 1753 } else { 1754 fits_copy_hdu (inputfptr, outfptr, 0, &stat); 1755 fits_copy_hdu (inputfptr, outfptr2, 0, &stat); 1756 } 1757 1758 fits_movrel_hdu (inputfptr, 1, NULL, &stat); 1759 extnum++; 1760 } 1761 1762 1763 if (stat == END_OF_FILE) stat = 0; 1764 1765 fits_close_file (outfptr2, &stat); 1766 fits_close_file (outfptr, &stat); 1767 fits_close_file (inputfptr, &stat); 1768 1769 if (stat) { 1770 fits_report_error (stderr, stat); 1771 } 1772 return(0); 1773 } 1774 /*--------------------------------------------------------------------------*/ 1775 int fp_pack_hdu (fitsfile *infptr, fitsfile *outfptr, fpstate fpvar, 1776 int *islossless, int *status) 1777 { 1778 fitsfile *tempfile; 1779 long naxes[9] = {1, 1, 1, 1, 1, 1, 1, 1, 1}; 1780 int stat=0, totpix=0, naxis=0, ii, hdutype, bitpix; 1781 int tstatus, hdunum; 1782 double bscale, rescale; 1783 1784 char outfits[SZ_STR], fzalgor[FLEN_VALUE]; 1785 long headstart, datastart, dataend, datasize; 1786 double noisemin; 1787 /* structure to hold image statistics (defined in fpack.h) */ 1788 imgstats imagestats; 1789 1790 if (*status) return(0); 1791 1792 fits_get_hdu_type (infptr, &hdutype, &stat); 1793 1794 if (hdutype == IMAGE_HDU) { 1795 fits_get_img_param (infptr, 9, &bitpix, &naxis, naxes, &stat); 1796 for (totpix=1, ii=0; ii < 9; ii++) totpix *= naxes[ii]; 1797 } 1798 1799 /* check directive keyword to see if this HDU should not be compressed */ 1800 tstatus = 0; 1801 if (!fits_read_key(infptr, TSTRING, "FZALGOR", fzalgor, NULL, &tstatus) ) { 1802 if (!strcmp(fzalgor, "NONE") || !strcmp(fzalgor, "none") ) { 1803 fits_copy_hdu (infptr, outfptr, 0, &stat); 1804 1805 *status = stat; 1806 return(0); 1807 } 1808 } 1809 1810 /* =============================================================== */ 1811 /* This block is only for binary table compression */ 1812 if (hdutype == BINARY_TBL && fpvar.do_tables) { 1813 1814 fits_get_hduaddr(infptr, &headstart, &datastart, &dataend, status); 1815 datasize = dataend - datastart; 1816 1817 if (datasize <= 2880) { 1818 /* data is less than 1 FITS block in size, so don't compress */ 1819 fits_copy_hdu (infptr, outfptr, 0, &stat); 1820 } else { 1821 fits_compress_table (infptr, outfptr, &stat); 1822 } 1823 1824 *status = stat; 1825 return(0); 1826 } 1827 /* =============================================================== */ 1828 1829 /* If this is not a non-null image HDU, just copy it verbatim */ 1830 if (fits_is_compressed_image (infptr, &stat) || hdutype != IMAGE_HDU || 1831 naxis == 0 || totpix == 0 || !fpvar.do_images) { 1832 fits_copy_hdu (infptr, outfptr, 0, &stat); 1833 1834 } else { /* remaining code deals only with IMAGE HDUs */ 1835 1836 /* special case: rescale a scaled integer image to reduce noise? */ 1837 if (fpvar.rescale_noise != 0. && bitpix > 0 && bitpix < LONGLONG_IMG) { 1838 1839 tstatus = 0; 1840 fits_read_key(infptr, TDOUBLE, "BSCALE", &bscale, 0, &tstatus); 1841 if (tstatus == 0 && bscale != 1.0) { /* image must be scaled */ 1842 1843 if (bitpix == LONG_IMG) 1844 fp_i4stat(infptr, naxis, naxes, &imagestats, &stat); 1845 else 1846 fp_i2stat(infptr, naxis, naxes, &imagestats, &stat); 1847 1848 /* use the minimum of the MAD 2nd, 3rd, and 5th order noise estimates */ 1849 noisemin = imagestats.noise3; 1850 if (imagestats.noise2 != 0. && imagestats.noise2 < noisemin) noisemin = imagestats.noise2; 1851 if (imagestats.noise5 != 0. && imagestats.noise5 < noisemin) noisemin = imagestats.noise5; 1852 1853 rescale = noisemin / fpvar.rescale_noise; 1854 if (rescale > 1.0) { 1855 1856 /* all the criteria are met, so create a temporary file that */ 1857 /* contains a rescaled version of the image, in output directory */ 1858 1859 /* create temporary file name */ 1860 fits_file_name(outfptr, outfits, &stat); /* get the output file name */ 1861 fp_tmpnam("Tmp3", outfits, tempfilename3); 1862 1863 fits_create_file(&tempfile, tempfilename3, &stat); 1864 1865 fits_get_hdu_num(infptr, &hdunum); 1866 if (hdunum != 1) { 1867 1868 /* the input hdu is an image extension, so create dummy primary */ 1869 fits_create_img(tempfile, 8, 0, naxes, &stat); 1870 } 1871 1872 fits_copy_header(infptr, tempfile, &stat); /* copy the header */ 1873 1874 /* rescale the data, so that it will compress more efficiently */ 1875 if (bitpix == LONG_IMG) 1876 fp_i4rescale(infptr, naxis, naxes, rescale, tempfile, &stat); 1877 else 1878 fp_i2rescale(infptr, naxis, naxes, rescale, tempfile, &stat); 1879 1880 1881 /* scale the BSCALE keyword by the inverse factor */ 1882 1883 bscale = bscale * rescale; 1884 fits_update_key(tempfile, TDOUBLE, "BSCALE", &bscale, 0, &stat); 1885 1886 /* rescan the header, to reset the actual scaling parameters */ 1887 fits_set_hdustruc(tempfile, &stat); 1888 1889 fits_img_compress (tempfile, outfptr, &stat); 1890 fits_delete_file (tempfile, &stat); 1891 tempfilename3[0] = '\0'; /* clear the temp filename */ 1892 *islossless = 0; /* used a lossy compression method */ 1893 1894 *status = stat; 1895 return(0); 1896 } 1897 } 1898 } 1899 1900 /* if requested to do lossy compression of integer images (by */ 1901 /* converting to float), then check if this HDU qualifies */ 1902 if ( (bitpix > 0) && (fpvar.int_to_float != 0) ) { 1903 1904 if (bitpix >= LONG_IMG) 1905 fp_i4stat(infptr, naxis, naxes, &imagestats, &stat); 1906 else 1907 fp_i2stat(infptr, naxis, naxes, &imagestats, &stat); 1908 1909 /* rescan the image header to reset scaling values (changed by fp_iNstat) */ 1910 ffrhdu(infptr, &hdutype, &stat); 1911 1912 /* use the minimum of the MAD 2nd, 3rd, and 5th order noise estimates */ 1913 noisemin = imagestats.noise3; 1914 if (imagestats.noise2 != 0. && imagestats.noise2 < noisemin) noisemin = imagestats.noise2; 1915 if (imagestats.noise5 != 0. && imagestats.noise5 < noisemin) noisemin = imagestats.noise5; 1916 1917 if ( (noisemin < (fpvar.n3ratio * fpvar.quantize_level) ) || 1918 (imagestats.noise3 < fpvar.n3min)) { 1919 1920 /* image contains too little noise to quantize effectively */ 1921 fits_set_lossy_int (outfptr, 0, &stat); 1922 1923 fits_get_hdu_num(infptr, &hdunum); 1924 1925 printf(" HDU %d does not meet noise criteria to be quantized, so losslessly compressed.\n", hdunum); 1926 1927 } else { 1928 /* compressed image is not identical to original */ 1929 *islossless = 0; 1930 } 1931 } 1932 1933 /* finally, do the actual image compression */ 1934 fits_img_compress (infptr, outfptr, &stat); 1935 1936 if (bitpix < 0 || 1937 (fpvar.comptype == HCOMPRESS_1 && fpvar.scale != 0.)) { 1938 1939 /* compressed image is not identical to original */ 1940 *islossless = 0; 1941 } 1942 } 1943 1944 *status = stat; 1945 return(0); 1946 } 1947 1948 /*--------------------------------------------------------------------------*/ 1949 int fp_unpack_hdu (fitsfile *infptr, fitsfile *outfptr, fpstate fpvar, int *status) 1950 { 1951 UNUSED(fpvar); 1952 1953 int hdutype, lval; 1954 1955 if (*status > 0) return(0); 1956 1957 fits_get_hdu_type (infptr, &hdutype, status); 1958 1959 /* =============================================================== */ 1960 /* This block is only for beta testing of binary table compression */ 1961 if (hdutype == BINARY_TBL) { 1962 1963 fits_read_key(infptr, TLOGICAL, "ZTABLE", &lval, NULL, status); 1964 1965 if (*status == 0 && lval != 0) { 1966 /* uncompress the table */ 1967 fits_uncompress_table (infptr, outfptr, status); 1968 } else { 1969 if (*status == KEY_NO_EXIST) /* table is not compressed */ 1970 *status = 0; 1971 fits_copy_hdu (infptr, outfptr, 0, status); 1972 } 1973 1974 return(0); 1975 /* =============================================================== */ 1976 1977 } else if (fits_is_compressed_image (infptr, status)) { 1978 /* uncompress the compressed image HDU */ 1979 fits_img_decompress (infptr, outfptr, status); 1980 } else { 1981 /* not a compressed image HDU, so just copy it to the output */ 1982 fits_copy_hdu (infptr, outfptr, 0, status); 1983 } 1984 1985 return(0); 1986 } 1987 /*--------------------------------------------------------------------------*/ 1988 int fits_read_image_speed (fitsfile *infptr, float *whole_elapse, 1989 float *whole_cpu, float *row_elapse, float *row_cpu, int *status) 1990 { 1991 unsigned char *carray, cnull = 0; 1992 short *sarray, snull=0; 1993 int bitpix, naxis, anynull, *iarray, inull = 0; 1994 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}; 1995 long inc[9]={1,1,1,1,1,1,1,1,1} ; 1996 float *earray, enull = 0, filesize; 1997 double *darray, dnull = 0; 1998 1999 if (*status) return(*status); 2000 2001 fits_get_img_param (infptr, 9, &bitpix, &naxis, naxes, status); 2002 2003 if (naxis != 2)return(*status); 2004 2005 lpixel[0] = naxes[0]; 2006 lpixel[1] = naxes[1]; 2007 2008 /* filesize in MB */ 2009 filesize = (float) (naxes[0] * abs(bitpix) / 8000000. * naxes[1]); 2010 2011 /* measure time required to read the raw image */ 2012 fits_set_bscale(infptr, 1.0, 0.0, status); 2013 *whole_elapse = 0.; 2014 *whole_cpu = 0; 2015 2016 if (bitpix == BYTE_IMG) { 2017 carray = calloc(naxes[1]*naxes[0], sizeof(char)); 2018 2019 /* remove any cached uncompressed tile 2020 (dangerous to directly modify the structure!) */ 2021 /* (infptr->Fptr)->tilerow = 0; */ 2022 2023 marktime(status); 2024 fits_read_subset(infptr, TBYTE, fpixel, lpixel, inc, &cnull, 2025 carray, &anynull, status); 2026 2027 /* get elapsped times */ 2028 gettime(whole_elapse, whole_cpu, status); 2029 2030 /* now read the image again, row by row */ 2031 if (row_elapse) { 2032 2033 /* remove any cached uncompressed tile 2034 (dangerous to directly modify the structure!) */ 2035 /* (infptr->Fptr)->tilerow = 0; */ 2036 2037 marktime(status); 2038 for (ii = 0; ii < naxes[1]; ii++) { 2039 fpixel[1] = ii+1; 2040 fits_read_pix(infptr, TBYTE, fpixel, naxes[0], &cnull, 2041 carray, &anynull, status); 2042 } 2043 /* get elapsped times */ 2044 gettime(row_elapse, row_cpu, status); 2045 } 2046 free(carray); 2047 2048 } else if (bitpix == SHORT_IMG) { 2049 sarray = calloc(naxes[0]*naxes[1], sizeof(short)); 2050 2051 marktime(status); 2052 fits_read_subset(infptr, TSHORT, fpixel, lpixel, inc, &snull, 2053 sarray, &anynull, status); 2054 2055 gettime(whole_elapse, whole_cpu, status); /* get elapsped times */ 2056 2057 /* now read the image again, row by row */ 2058 if (row_elapse) { 2059 marktime(status); 2060 for (ii = 0; ii < naxes[1]; ii++) { 2061 2062 fpixel[1] = ii+1; 2063 fits_read_pix(infptr, TSHORT, fpixel, naxes[0], &snull, 2064 sarray, &anynull, status); 2065 } 2066 /* get elapsped times */ 2067 gettime(row_elapse, row_cpu, status); 2068 } 2069 2070 free(sarray); 2071 2072 } else if (bitpix == LONG_IMG) { 2073 iarray = calloc(naxes[0]*naxes[1], sizeof(int)); 2074 2075 marktime(status); 2076 2077 fits_read_subset(infptr, TINT, fpixel, lpixel, inc, &inull, 2078 iarray, &anynull, status); 2079 2080 /* get elapsped times */ 2081 gettime(whole_elapse, whole_cpu, status); 2082 2083 2084 /* now read the image again, row by row */ 2085 if (row_elapse) { 2086 marktime(status); 2087 for (ii = 0; ii < naxes[1]; ii++) { 2088 fpixel[1] = ii+1; 2089 fits_read_pix(infptr, TINT, fpixel, naxes[0], &inull, 2090 iarray, &anynull, status); 2091 } 2092 /* get elapsped times */ 2093 gettime(row_elapse, row_cpu, status); 2094 } 2095 2096 2097 free(iarray); 2098 2099 } else if (bitpix == FLOAT_IMG) { 2100 earray = calloc(naxes[1]*naxes[0], sizeof(float)); 2101 2102 marktime(status); 2103 2104 fits_read_subset(infptr, TFLOAT, fpixel, lpixel, inc, &enull, 2105 earray, &anynull, status); 2106 2107 /* get elapsped times */ 2108 gettime(whole_elapse, whole_cpu, status); 2109 2110 /* now read the image again, row by row */ 2111 if (row_elapse) { 2112 marktime(status); 2113 for (ii = 0; ii < naxes[1]; ii++) { 2114 fpixel[1] = ii+1; 2115 fits_read_pix(infptr, TFLOAT, fpixel, naxes[0], &enull, 2116 earray, &anynull, status); 2117 } 2118 /* get elapsped times */ 2119 gettime(row_elapse, row_cpu, status); 2120 } 2121 2122 free(earray); 2123 2124 } else if (bitpix == DOUBLE_IMG) { 2125 darray = calloc(naxes[1]*naxes[0], sizeof(double)); 2126 2127 marktime(status); 2128 2129 fits_read_subset(infptr, TDOUBLE, fpixel, lpixel, inc, &dnull, 2130 darray, &anynull, status); 2131 2132 /* get elapsped times */ 2133 gettime(whole_elapse, whole_cpu, status); 2134 2135 /* now read the image again, row by row */ 2136 if (row_elapse) { 2137 marktime(status); 2138 for (ii = 0; ii < naxes[1]; ii++) { 2139 fpixel[1] = ii+1; 2140 fits_read_pix(infptr, TDOUBLE, fpixel, naxes[0], &dnull, 2141 darray, &anynull, status); 2142 } 2143 /* get elapsped times */ 2144 gettime(row_elapse, row_cpu, status); 2145 } 2146 2147 free(darray); 2148 } 2149 2150 if (whole_elapse) *whole_elapse = *whole_elapse / filesize; 2151 if (row_elapse) *row_elapse = *row_elapse / filesize; 2152 if (whole_cpu) *whole_cpu = *whole_cpu / filesize; 2153 if (row_cpu) *row_cpu = *row_cpu / filesize; 2154 2155 return(*status); 2156 } 2157 /*--------------------------------------------------------------------------*/ 2158 int fp_test_hdu (fitsfile *infptr, fitsfile *outfptr, fitsfile *outfptr2, 2159 fpstate fpvar, int *status) 2160 { 2161 /* This routine is only used for performance testing of image HDUs. */ 2162 /* Use fp_test_table for testing binary table HDUs. */ 2163 2164 int stat = 0, hdutype, comptype; 2165 char ctype[20], lossless[4]; 2166 long headstart, datastart, dataend; 2167 float origdata = 0., compressdata = 0.; 2168 float compratio = 0., packcpu = 0., unpackcpu = 0.; 2169 float elapse, whole_elapse, row_elapse, whole_cpu, row_cpu; 2170 unsigned long datasum1, datasum2, hdusum; 2171 2172 if (*status) return(0); 2173 2174 origdata = 0; 2175 compressdata = 0; 2176 compratio = 0.; 2177 lossless[0] = '\0'; 2178 2179 fits_get_compression_type(outfptr, &comptype, &stat); 2180 if (comptype == RICE_1) 2181 strcpy(ctype, "RICE"); 2182 else if (comptype == GZIP_1) 2183 strcpy(ctype, "GZIP1"); 2184 else if (comptype == GZIP_2) 2185 strcpy(ctype, "GZIP2");/* 2186 else if (comptype == BZIP2_1) 2187 strcpy(ctype, "BZIP2"); 2188 */ 2189 else if (comptype == PLIO_1) 2190 strcpy(ctype, "PLIO"); 2191 else if (comptype == HCOMPRESS_1) 2192 strcpy(ctype, "HCOMP"); 2193 else if (comptype == NOCOMPRESS) 2194 strcpy(ctype, "NONE"); 2195 else { 2196 fp_msg ("Error: unsupported image compression type "); 2197 *status = DATA_COMPRESSION_ERR; 2198 return(0); 2199 } 2200 2201 /* -------------- COMPRESS the image ------------------ */ 2202 2203 marktime(&stat); 2204 2205 fits_img_compress (infptr, outfptr, &stat); 2206 2207 /* get elapsped times */ 2208 gettime(&elapse, &packcpu, &stat); 2209 2210 /* get elapsed and cpu times need to read the compressed image */ 2211 fits_read_image_speed (outfptr, &whole_elapse, &whole_cpu, 2212 &row_elapse, &row_cpu, &stat); 2213 2214 if (!stat) { 2215 2216 /* -------------- UNCOMPRESS the image ------------------ */ 2217 2218 /* remove any cached uncompressed tile 2219 (dangerous to directly modify the structure!) */ 2220 /* (outfptr->Fptr)->tilerow = 0; */ 2221 marktime(&stat); 2222 2223 fits_img_decompress (outfptr, outfptr2, &stat); 2224 2225 /* get elapsped times */ 2226 gettime(&elapse, &unpackcpu, &stat); 2227 2228 /* ----------------------------------------------------- */ 2229 2230 2231 /* get sizes of original and compressed images */ 2232 2233 fits_get_hduaddr(infptr, &headstart, &datastart, &dataend, &stat); 2234 origdata = (float) ((dataend - datastart)/1000000.); 2235 2236 fits_get_hduaddr(outfptr, &headstart, &datastart, &dataend, &stat); 2237 compressdata = (float) ((dataend - datastart)/1000000.); 2238 2239 if (compressdata != 0) 2240 compratio = (float) origdata / (float) compressdata; 2241 2242 /* is this uncompressed image identical to the original? */ 2243 2244 fits_get_chksum(infptr, &datasum1, &hdusum, &stat); 2245 fits_get_chksum(outfptr2, &datasum2, &hdusum, &stat); 2246 2247 if ( datasum1 == datasum2) { 2248 strcpy(lossless, "Yes"); 2249 } else { 2250 strcpy(lossless, "No"); 2251 } 2252 2253 printf(" %-5s %6.2f %7.2f ->%7.2f %7.2f %7.2f %s %5.3f %5.3f %5.3f %5.3f\n", 2254 ctype, compratio, origdata, compressdata, 2255 packcpu, unpackcpu, lossless, whole_elapse, whole_cpu, 2256 row_elapse, row_cpu); 2257 2258 2259 if (fpvar.outfile[0]) { 2260 fprintf(outreport," %6.3f %5.2f %5.2f %s %7.3f %7.3f %7.3f %7.3f", 2261 compratio, packcpu, unpackcpu, lossless, whole_elapse, whole_cpu, 2262 row_elapse, row_cpu); 2263 } 2264 2265 /* delete the output HDUs to concerve disk space */ 2266 2267 fits_delete_hdu(outfptr, &hdutype, &stat); 2268 fits_delete_hdu(outfptr2, &hdutype, &stat); 2269 2270 } else { 2271 printf(" %-5s (unable to compress image)\n", ctype); 2272 } 2273 2274 /* try to recover from any compression errors */ 2275 if (stat == DATA_COMPRESSION_ERR) stat = 0; 2276 2277 *status = stat; 2278 return(0); 2279 } 2280 /*--------------------------------------------------------------------------*/ 2281 int fp_test_table (fitsfile *infptr, fitsfile *outfptr, fitsfile *outfptr2, 2282 fpstate fpvar, int *status) 2283 { 2284 /* this routine is for performance testing of the table compression methods */ 2285 UNUSED(fpvar); 2286 UNUSED(outfptr2); 2287 2288 int stat = 0, hdutype, tstatus = 0; 2289 char fzalgor[FLEN_VALUE]; 2290 LONGLONG headstart, datastart, dataend; 2291 float elapse, cpu; 2292 2293 if (*status) return(0); 2294 2295 /* check directive keyword to see if this HDU should not be compressed */ 2296 if (!fits_read_key(infptr, TSTRING, "FZALGOR", fzalgor, NULL, &tstatus) ) { 2297 if (!strcmp(fzalgor, "NONE") || !strcmp(fzalgor, "none")) { 2298 return(0); 2299 } 2300 } 2301 2302 fits_get_hduaddrll(infptr, &headstart, &datastart, &dataend, status); 2303 2304 /* can't compress small tables with less than 2880 bytes of data */ 2305 if (dataend - datastart <= 2880) { 2306 return(0); 2307 } 2308 2309 marktime(&stat); 2310 stat= -999; /* set special flag value */ 2311 fits_compress_table (infptr, outfptr, &stat); 2312 2313 /* get elapsped times */ 2314 gettime(&elapse, &cpu, &stat); 2315 2316 fits_delete_hdu(outfptr, &hdutype, &stat); 2317 2318 printf("\nElapsed time = %f, cpu = %f\n", elapse, cpu); 2319 2320 fits_report_error (stderr, stat); 2321 2322 return(0); 2323 } 2324 /*--------------------------------------------------------------------------*/ 2325 int marktime(int *status) 2326 { 2327 #if defined(unix) || defined(__unix__) || defined(__unix) 2328 struct timeval tv; 2329 /* struct timezone tz; */ 2330 2331 /* gettimeofday (&tv, &tz); */ 2332 gettimeofday (&tv, NULL); 2333 2334 startsec = tv.tv_sec; 2335 startmilli = tv.tv_usec/1000; 2336 2337 scpu = clock(); 2338 #else 2339 /* don't support high timing precision on Windows machines */ 2340 startsec = 0; 2341 startmilli = 0; 2342 2343 scpu = clock(); 2344 #endif 2345 return( *status ); 2346 } 2347 /*--------------------------------------------------------------------------*/ 2348 int gettime(float *elapse, float *elapscpu, int *status) 2349 { 2350 #if defined(unix) || defined(__unix__) || defined(__unix) 2351 struct timeval tv; 2352 /* struct timezone tz; */ 2353 int stopmilli; 2354 long stopsec; 2355 2356 /* gettimeofday (&tv, &tz); */ 2357 gettimeofday (&tv, NULL); 2358 ecpu = clock(); 2359 2360 stopmilli = tv.tv_usec/1000; 2361 stopsec = tv.tv_sec; 2362 2363 *elapse = (stopsec - startsec) + (stopmilli - startmilli)/1000.; 2364 *elapscpu = (ecpu - scpu) * 1.0 / CLOCKTICKS; 2365 /* 2366 printf(" (start: %ld + %d), stop: (%ld + %d) elapse: %f\n ", 2367 startsec,startmilli,stopsec, stopmilli, *elapse); 2368 */ 2369 #else 2370 /* set the elapsed time the same as the CPU time on Windows machines */ 2371 *elapscpu = (float) ((ecpu - scpu) * 1.0 / CLOCKTICKS); 2372 *elapse = *elapscpu; 2373 #endif 2374 return( *status ); 2375 } 2376 /*--------------------------------------------------------------------------*/ 2377 int fp_i2stat(fitsfile *infptr, int naxis, long *naxes, imgstats *imagestats, int *status) 2378 { 2379 /* 2380 read the central XSAMPLE by YSAMPLE region of pixels in the int*2 image, 2381 and then compute basic statistics: min, max, mean, sigma, mean diff, etc. 2382 */ 2383 2384 long fpixel[9] = {1,1,1,1,1,1,1,1,1}; 2385 long lpixel[9] = {1,1,1,1,1,1,1,1,1}; 2386 long inc[9] = {1,1,1,1,1,1,1,1,1}; 2387 long i1, i2, npix, ngood, nx, ny; 2388 short *intarray, minvalue, maxvalue, nullvalue; 2389 int anynul, tstatus, checknull = 1; 2390 double mean, sigma, noise1, noise2, noise3, noise5; 2391 2392 /* select the middle XSAMPLE by YSAMPLE area of the image */ 2393 i1 = naxes[0]/2 - (XSAMPLE/2 - 1); 2394 i2 = naxes[0]/2 + (XSAMPLE/2); 2395 if (i1 < 1) i1 = 1; 2396 if (i2 > naxes[0]) i2 = naxes[0]; 2397 fpixel[0] = i1; 2398 lpixel[0] = i2; 2399 nx = i2 - i1 +1; 2400 2401 if (naxis > 1) { 2402 i1 = naxes[1]/2 - (YSAMPLE/2 - 1); 2403 i2 = naxes[1]/2 + (YSAMPLE/2); 2404 if (i1 < 1) i1 = 1; 2405 if (i2 > naxes[1]) i2 = naxes[1]; 2406 fpixel[1] = i1; 2407 lpixel[1] = i2; 2408 } 2409 ny = i2 - i1 +1; 2410 2411 npix = nx * ny; 2412 2413 /* if there are higher dimensions, read the middle plane of the cube */ 2414 if (naxis > 2) { 2415 fpixel[2] = naxes[2]/2 + 1; 2416 lpixel[2] = naxes[2]/2 + 1; 2417 } 2418 2419 intarray = calloc(npix, sizeof(short)); 2420 if (!intarray) { 2421 *status = MEMORY_ALLOCATION; 2422 return(*status); 2423 } 2424 2425 /* turn off any scaling of the integer pixel values */ 2426 fits_set_bscale(infptr, 1.0, 0.0, status); 2427 2428 fits_read_subset_sht(infptr, 0, naxis, naxes, fpixel, lpixel, inc, 2429 0, intarray, &anynul, status); 2430 2431 /* read the null value keyword (BLANK) if present */ 2432 tstatus = 0; 2433 fits_read_key(infptr, TSHORT, "BLANK", &nullvalue, 0, &tstatus); 2434 if (tstatus) { 2435 nullvalue = 0; 2436 checknull = 0; 2437 } 2438 2439 /* compute statistics of the image */ 2440 2441 fits_img_stats_short(intarray, nx, ny, checknull, nullvalue, 2442 &ngood, &minvalue, &maxvalue, &mean, &sigma, &noise1, &noise2, &noise3, &noise5, status); 2443 2444 imagestats->n_nulls = npix - ngood; 2445 imagestats->minval = minvalue; 2446 imagestats->maxval = maxvalue; 2447 imagestats->mean = mean; 2448 imagestats->sigma = sigma; 2449 imagestats->noise1 = noise1; 2450 imagestats->noise2 = noise2; 2451 imagestats->noise3 = noise3; 2452 imagestats->noise5 = noise5; 2453 2454 free(intarray); 2455 return(*status); 2456 } 2457 /*--------------------------------------------------------------------------*/ 2458 int fp_i4stat(fitsfile *infptr, int naxis, long *naxes, imgstats *imagestats, int *status) 2459 { 2460 /* 2461 read the central XSAMPLE by YSAMPLE region of pixels in the int*2 image, 2462 and then compute basic statistics: min, max, mean, sigma, mean diff, etc. 2463 */ 2464 2465 long fpixel[9] = {1,1,1,1,1,1,1,1,1}; 2466 long lpixel[9] = {1,1,1,1,1,1,1,1,1}; 2467 long inc[9] = {1,1,1,1,1,1,1,1,1}; 2468 long i1, i2, npix, ngood, nx, ny; 2469 int *intarray, minvalue, maxvalue, nullvalue; 2470 int anynul, tstatus, checknull = 1; 2471 double mean, sigma, noise1, noise2, noise3, noise5; 2472 2473 /* select the middle XSAMPLE by YSAMPLE area of the image */ 2474 i1 = naxes[0]/2 - (XSAMPLE/2 - 1); 2475 i2 = naxes[0]/2 + (XSAMPLE/2); 2476 if (i1 < 1) i1 = 1; 2477 if (i2 > naxes[0]) i2 = naxes[0]; 2478 fpixel[0] = i1; 2479 lpixel[0] = i2; 2480 nx = i2 - i1 +1; 2481 2482 if (naxis > 1) { 2483 i1 = naxes[1]/2 - (YSAMPLE/2 - 1); 2484 i2 = naxes[1]/2 + (YSAMPLE/2); 2485 if (i1 < 1) i1 = 1; 2486 if (i2 > naxes[1]) i2 = naxes[1]; 2487 fpixel[1] = i1; 2488 lpixel[1] = i2; 2489 } 2490 ny = i2 - i1 +1; 2491 2492 npix = nx * ny; 2493 2494 /* if there are higher dimensions, read the middle plane of the cube */ 2495 if (naxis > 2) { 2496 fpixel[2] = naxes[2]/2 + 1; 2497 lpixel[2] = naxes[2]/2 + 1; 2498 } 2499 2500 intarray = calloc(npix, sizeof(int)); 2501 if (!intarray) { 2502 *status = MEMORY_ALLOCATION; 2503 return(*status); 2504 } 2505 2506 /* turn off any scaling of the integer pixel values */ 2507 fits_set_bscale(infptr, 1.0, 0.0, status); 2508 2509 fits_read_subset_int(infptr, 0, naxis, naxes, fpixel, lpixel, inc, 2510 0, intarray, &anynul, status); 2511 2512 /* read the null value keyword (BLANK) if present */ 2513 tstatus = 0; 2514 fits_read_key(infptr, TINT, "BLANK", &nullvalue, 0, &tstatus); 2515 if (tstatus) { 2516 nullvalue = 0; 2517 checknull = 0; 2518 } 2519 2520 /* compute statistics of the image */ 2521 2522 fits_img_stats_int(intarray, nx, ny, checknull, nullvalue, 2523 &ngood, &minvalue, &maxvalue, &mean, &sigma, &noise1, &noise2, &noise3, &noise5, status); 2524 2525 imagestats->n_nulls = npix - ngood; 2526 imagestats->minval = minvalue; 2527 imagestats->maxval = maxvalue; 2528 imagestats->mean = mean; 2529 imagestats->sigma = sigma; 2530 imagestats->noise1 = noise1; 2531 imagestats->noise2 = noise2; 2532 imagestats->noise3 = noise3; 2533 imagestats->noise5 = noise5; 2534 2535 free(intarray); 2536 return(*status); 2537 } 2538 /*--------------------------------------------------------------------------*/ 2539 int fp_r4stat(fitsfile *infptr, int naxis, long *naxes, imgstats *imagestats, int *status) 2540 { 2541 /* 2542 read the central XSAMPLE by YSAMPLE region of pixels in the int*2 image, 2543 and then compute basic statistics: min, max, mean, sigma, mean diff, etc. 2544 */ 2545 2546 long fpixel[9] = {1,1,1,1,1,1,1,1,1}; 2547 long lpixel[9] = {1,1,1,1,1,1,1,1,1}; 2548 long inc[9] = {1,1,1,1,1,1,1,1,1}; 2549 long i1, i2, npix, ngood, nx, ny; 2550 float *array, minvalue, maxvalue, nullvalue = FLOATNULLVALUE; 2551 int anynul,checknull = 1; 2552 double mean, sigma, noise1, noise2, noise3, noise5; 2553 2554 /* select the middle XSAMPLE by YSAMPLE area of the image */ 2555 i1 = naxes[0]/2 - (XSAMPLE/2 - 1); 2556 i2 = naxes[0]/2 + (XSAMPLE/2); 2557 if (i1 < 1) i1 = 1; 2558 if (i2 > naxes[0]) i2 = naxes[0]; 2559 fpixel[0] = i1; 2560 lpixel[0] = i2; 2561 nx = i2 - i1 +1; 2562 2563 if (naxis > 1) { 2564 i1 = naxes[1]/2 - (YSAMPLE/2 - 1); 2565 i2 = naxes[1]/2 + (YSAMPLE/2); 2566 if (i1 < 1) i1 = 1; 2567 if (i2 > naxes[1]) i2 = naxes[1]; 2568 fpixel[1] = i1; 2569 lpixel[1] = i2; 2570 } 2571 ny = i2 - i1 +1; 2572 2573 npix = nx * ny; 2574 2575 /* if there are higher dimensions, read the middle plane of the cube */ 2576 if (naxis > 2) { 2577 fpixel[2] = naxes[2]/2 + 1; 2578 lpixel[2] = naxes[2]/2 + 1; 2579 } 2580 2581 array = calloc(npix, sizeof(float)); 2582 if (!array) { 2583 *status = MEMORY_ALLOCATION; 2584 return(*status); 2585 } 2586 2587 fits_read_subset_flt(infptr, 0, naxis, naxes, fpixel, lpixel, inc, 2588 nullvalue, array, &anynul, status); 2589 2590 /* are there any null values in the array? */ 2591 if (!anynul) { 2592 nullvalue = 0.; 2593 checknull = 0; 2594 } 2595 2596 /* compute statistics of the image */ 2597 2598 fits_img_stats_float(array, nx, ny, checknull, nullvalue, 2599 &ngood, &minvalue, &maxvalue, &mean, &sigma, &noise1, &noise2, &noise3, &noise5, status); 2600 2601 imagestats->n_nulls = npix - ngood; 2602 imagestats->minval = minvalue; 2603 imagestats->maxval = maxvalue; 2604 imagestats->mean = mean; 2605 imagestats->sigma = sigma; 2606 imagestats->noise1 = noise1; 2607 imagestats->noise2 = noise2; 2608 imagestats->noise3 = noise3; 2609 imagestats->noise5 = noise5; 2610 2611 free(array); 2612 return(*status); 2613 } 2614 /*--------------------------------------------------------------------------*/ 2615 int fp_i2rescale(fitsfile *infptr, int naxis, long *naxes, double rescale, 2616 fitsfile *outfptr, int *status) 2617 { 2618 /* 2619 divide the integer pixel values in the input file by rescale, 2620 and write back out to the output file.. 2621 */ 2622 2623 long ii, jj, nelem = 1, nx, ny; 2624 short *intarray, nullvalue; 2625 int anynul, tstatus, checknull = 1; 2626 2627 nx = naxes[0]; 2628 ny = 1; 2629 2630 for (ii = 1; ii < naxis; ii++) { 2631 ny = ny * naxes[ii]; 2632 } 2633 2634 intarray = calloc(nx, sizeof(short)); 2635 if (!intarray) { 2636 *status = MEMORY_ALLOCATION; 2637 return(*status); 2638 } 2639 2640 /* read the null value keyword (BLANK) if present */ 2641 tstatus = 0; 2642 fits_read_key(infptr, TSHORT, "BLANK", &nullvalue, 0, &tstatus); 2643 if (tstatus) { 2644 checknull = 0; 2645 } 2646 2647 /* turn off any scaling of the integer pixel values */ 2648 fits_set_bscale(infptr, 1.0, 0.0, status); 2649 fits_set_bscale(outfptr, 1.0, 0.0, status); 2650 2651 for (ii = 0; ii < ny; ii++) { 2652 2653 fits_read_img_sht(infptr, 1, nelem, nx, 2654 0, intarray, &anynul, status); 2655 2656 if (checknull) { 2657 for (jj = 0; jj < nx; jj++) { 2658 if (intarray[jj] != nullvalue) 2659 intarray[jj] = NSHRT( (intarray[jj] / rescale) ); 2660 } 2661 } else { 2662 for (jj = 0; jj < nx; jj++) 2663 intarray[jj] = NSHRT( (intarray[jj] / rescale) ); 2664 } 2665 2666 fits_write_img_sht(outfptr, 1, nelem, nx, intarray, status); 2667 2668 nelem += nx; 2669 } 2670 2671 free(intarray); 2672 return(*status); 2673 } 2674 /*--------------------------------------------------------------------------*/ 2675 int fp_i4rescale(fitsfile *infptr, int naxis, long *naxes, double rescale, 2676 fitsfile *outfptr, int *status) 2677 { 2678 /* 2679 divide the integer pixel values in the input file by rescale, 2680 and write back out to the output file.. 2681 */ 2682 2683 long ii, jj, nelem = 1, nx, ny; 2684 int *intarray, nullvalue; 2685 int anynul, tstatus, checknull = 1; 2686 2687 nx = naxes[0]; 2688 ny = 1; 2689 2690 for (ii = 1; ii < naxis; ii++) { 2691 ny = ny * naxes[ii]; 2692 } 2693 2694 intarray = calloc(nx, sizeof(int)); 2695 if (!intarray) { 2696 *status = MEMORY_ALLOCATION; 2697 return(*status); 2698 } 2699 2700 /* read the null value keyword (BLANK) if present */ 2701 tstatus = 0; 2702 fits_read_key(infptr, TINT, "BLANK", &nullvalue, 0, &tstatus); 2703 if (tstatus) { 2704 checknull = 0; 2705 } 2706 2707 /* turn off any scaling of the integer pixel values */ 2708 fits_set_bscale(infptr, 1.0, 0.0, status); 2709 fits_set_bscale(outfptr, 1.0, 0.0, status); 2710 2711 for (ii = 0; ii < ny; ii++) { 2712 2713 fits_read_img_int(infptr, 1, nelem, nx, 2714 0, intarray, &anynul, status); 2715 2716 if (checknull) { 2717 for (jj = 0; jj < nx; jj++) { 2718 if (intarray[jj] != nullvalue) 2719 intarray[jj] = NINT( (intarray[jj] / rescale) ); 2720 } 2721 } else { 2722 for (jj = 0; jj < nx; jj++) 2723 intarray[jj] = NINT( (intarray[jj] / rescale) ); 2724 } 2725 2726 fits_write_img_int(outfptr, 1, nelem, nx, intarray, status); 2727 2728 nelem += nx; 2729 } 2730 2731 free(intarray); 2732 return(*status); 2733 } 2734 /* ======================================================================== 2735 * Signal and error handler. 2736 */ 2737 void abort_fpack(int sig) 2738 { 2739 /* clean up by deleting temporary files */ 2740 UNUSED(sig); 2741 2742 if (tempfilename[0]) { 2743 remove(tempfilename); 2744 } 2745 if (tempfilename2[0]) { 2746 remove(tempfilename2); 2747 } 2748 if (tempfilename3[0]) { 2749 remove(tempfilename3); 2750 } 2751 exit(-1); 2752 }