LSSTApplications  17.0+124,17.0+14,17.0+73,18.0.0+37,18.0.0+80,18.0.0-4-g68ffd23+4,18.1.0-1-g0001055+12,18.1.0-1-g03d53ef+5,18.1.0-1-g1349e88+55,18.1.0-1-g2505f39+44,18.1.0-1-g5315e5e+4,18.1.0-1-g5e4b7ea+14,18.1.0-1-g7e8fceb+4,18.1.0-1-g85f8cd4+48,18.1.0-1-g8ff0b9f+4,18.1.0-1-ga2c679d+1,18.1.0-1-gd55f500+35,18.1.0-10-gb58edde+2,18.1.0-11-g0997b02+4,18.1.0-13-gfe4edf0b+12,18.1.0-14-g259bd21+21,18.1.0-19-gdb69f3f+2,18.1.0-2-g5f9922c+24,18.1.0-2-gd3b74e5+11,18.1.0-2-gfbf3545+32,18.1.0-26-g728bddb4+5,18.1.0-27-g6ff7ca9+2,18.1.0-3-g52aa583+25,18.1.0-3-g8ea57af+9,18.1.0-3-gb69f684+42,18.1.0-3-gfcaddf3+6,18.1.0-32-gd8786685a,18.1.0-4-gf3f9b77+6,18.1.0-5-g1dd662b+2,18.1.0-5-g6dbcb01+41,18.1.0-6-gae77429+3,18.1.0-7-g9d75d83+9,18.1.0-7-gae09a6d+30,18.1.0-9-gc381ef5+4,w.2019.45
LSSTDataManagementBasePackage
fits.cc
Go to the documentation of this file.
1 // -*- lsst-c++ -*-
2 
3 #include <cstdint>
4 #include <cstdio>
5 #include <complex>
6 #include <cmath>
7 #include <sstream>
8 #include <unordered_set>
9 #include <unordered_map>
10 
11 #include "fitsio.h"
12 extern "C" {
13 #include "fitsio2.h"
14 }
15 
16 #include "boost/regex.hpp"
17 #include "boost/filesystem.hpp"
18 #include "boost/preprocessor/seq/for_each.hpp"
19 #include "boost/format.hpp"
20 
21 #include "lsst/pex/exceptions.h"
22 #include "lsst/log/Log.h"
23 #include "lsst/afw/fits.h"
24 #include "lsst/geom/Angle.h"
25 #include "lsst/afw/geom/wcsUtils.h"
27 
28 namespace lsst {
29 namespace afw {
30 namespace fits {
31 
32 // ----------------------------------------------------------------------------------------------------------
33 // ---- Miscellaneous utilities -----------------------------------------------------------------------------
34 // ----------------------------------------------------------------------------------------------------------
35 
36 namespace {
37 
38 /*
39  * Format a PropertySet into a FITS header string using simplifying assumptions.
40  *
41  * See @ref makeLimitedFitsHeader for details.
42  *
43  * @param[in] paramNames Names of properties to format
44  * @param[in] metadata Metadata to format
45  * @return a FITS header string (exactly 80 characters per entry, no line terminators)
46  */
47 std::string makeLimitedFitsHeaderImpl(std::vector<std::string> const &paramNames,
48  daf::base::PropertySet const &metadata) {
50  for (auto const &fullName : paramNames) {
51  std::size_t lastPeriod = fullName.rfind(char('.'));
52  auto name = (lastPeriod == std::string::npos) ? fullName : fullName.substr(lastPeriod + 1);
53  std::type_info const &type = metadata.typeOf(name);
54 
55  std::string out = "";
56  out.reserve(80);
57  if (name.size() > 8) {
58  continue; // The name is too long for a FITS keyword; skip this item
59  }
60  out = (boost::format("%-8s= ") % name).str();
61 
62  if (type == typeid(bool)) {
63  out += metadata.get<bool>(name) ? "T" : "F";
64  } else if (type == typeid(std::uint8_t)) {
65  out += (boost::format("%20d") % static_cast<int>(metadata.get<std::uint8_t>(name))).str();
66  } else if (type == typeid(int)) {
67  out += (boost::format("%20d") % metadata.get<int>(name)).str();
68  } else if (type == typeid(double)) {
69  double value = metadata.get<double>(name);
70  if (!std::isnan(value)) {
71  // use G because FITS wants uppercase E for exponents
72  out += (boost::format("%#20.17G") % value).str();
73  } else {
74  LOGLS_WARN("afw.fits",
75  boost::format("In %s, found NaN in metadata item '%s'") %
76  BOOST_CURRENT_FUNCTION % name);
77  // Convert it to FITS undefined
78  out += " ";
79  }
80  } else if (type == typeid(float)) {
81  float value = metadata.get<float>(name);
82  if (!std::isnan(value)) {
83  out += (boost::format("%#20.15G") % value).str();
84  } else {
85  LOGLS_WARN("afw.fits",
86  boost::format("In %s, found NaN in metadata item '%s'") %
87  BOOST_CURRENT_FUNCTION % name);
88  // Convert it to FITS undefined
89  out += " ";
90  }
91  } else if (type == typeid(std::nullptr_t)) {
92  out += " ";
93  } else if (type == typeid(std::string)) {
94  out += "'" + metadata.get<std::string>(name) + "'";
95  if (out.size() > 80) {
96  continue; // Formatted data is too long; skip this item
97  }
98  }
99 
100  int const len = out.size();
101  if (len < 80) {
102  out += std::string(80 - len, ' ');
103  } else if (len > 80) {
104  // non-string item has a formatted value that is too long; this should never happen
105  throw LSST_EXCEPT(pex::exceptions::LogicError,
106  "Formatted data too long: " + std::to_string(len) + " > 80: \"" + out + "\"");
107  }
108 
109  result << out;
110  }
111 
112  return result.str();
113 }
114 
123 class StringStartSet {
124 public:
126  StringStartSet(std::initializer_list<std::string> const &input) : _minSize(-1) {
127  for (auto const &word : input) {
128  std::size_t const size = word.size();
129  if (size < _minSize) {
130  _minSize = size;
131  }
132  }
133  for (auto const &word : input) {
134  std::string const start = startString(word);
135  assert(_words.count(start) == 0); // This should be the only word that starts this way
136  _words[start] = word;
137  }
138  }
139 
141  bool matches(std::string const &key) const {
142  auto const iter = _words.find(startString(key));
143  if (iter == _words.end()) {
144  return false;
145  }
146  // Check that the full word matches too
147  std::string const &word = iter->second;
148  return key.compare(0, word.size(), word) == 0;
149  }
150 
151 private:
153 
155  std::string startString(std::string const &word) const { return word.substr(0, _minSize); }
156 
157  std::size_t _minSize; // Minimum length of provided words
158  Map _words; // Start of words --> full word
159 };
160 
166 static std::unordered_set<std::string> const ignoreKeys = {
167  // FITS core keywords
168  "SIMPLE", "BITPIX", "NAXIS", "EXTEND", "GCOUNT", "PCOUNT", "XTENSION", "TFIELDS", "BSCALE", "BZERO",
169  // FITS compression keywords
170  "ZBITPIX", "ZIMAGE", "ZCMPTYPE", "ZSIMPLE", "ZEXTEND", "ZBLANK", "ZDATASUM", "ZHECKSUM", "ZQUANTIZ",
171  // Not essential, but will prevent fitsverify warnings
172  "DATASUM", "CHECKSUM"};
173 
179 StringStartSet const ignoreKeyStarts{// FITS core keywords
180  "NAXIS", "TZERO", "TSCAL",
181  // FITS compression keywords
182  "ZNAXIS", "ZTILE", "ZNAME", "ZVAL"};
183 
189 StringStartSet const ignoreKeyStartsWrite{"TFORM", "TTYPE"};
190 
191 // Strip leading and trailing single quotes and whitespace from a string.
192 std::string strip(std::string const &s) {
193  if (s.empty()) return s;
194  std::size_t i1 = s.find_first_not_of(" '");
195  std::size_t i2 = s.find_last_not_of(" '");
196  return s.substr(i1, (i1 == std::string::npos) ? 0 : 1 + i2 - i1);
197 }
198 
199 // ---- FITS binary table format codes for various C++ types. -----------------------------------------------
200 
201 char getFormatCode(bool *) { return 'X'; }
202 char getFormatCode(std::string *) { return 'A'; }
203 char getFormatCode(std::int8_t *) { return 'S'; }
204 char getFormatCode(std::uint8_t *) { return 'B'; }
205 char getFormatCode(std::int16_t *) { return 'I'; }
206 char getFormatCode(std::uint16_t *) { return 'U'; }
207 char getFormatCode(std::int32_t *) { return 'J'; }
208 char getFormatCode(std::uint32_t *) { return 'V'; }
209 char getFormatCode(std::int64_t *) { return 'K'; }
210 char getFormatCode(float *) { return 'E'; }
211 char getFormatCode(double *) { return 'D'; }
212 char getFormatCode(std::complex<float> *) { return 'C'; }
213 char getFormatCode(std::complex<double> *) { return 'M'; }
214 char getFormatCode(lsst::geom::Angle *) { return 'D'; }
215 
216 // ---- Create a TFORM value for the given type and size ----------------------------------------------------
217 
218 template <typename T>
219 std::string makeColumnFormat(int size = 1) {
220  if (size > 0) {
221  return (boost::format("%d%c") % size % getFormatCode((T *)0)).str();
222  } else if (size < 0) {
223  // variable length, max size given as -size
224  return (boost::format("1Q%c(%d)") % getFormatCode((T *)0) % (-size)).str();
225  } else {
226  // variable length, max size unknown
227  return (boost::format("1Q%c") % getFormatCode((T *)0)).str();
228  }
229 }
230 
231 // ---- Traits class to get cfitsio type constants from templates -------------------------------------------
232 
233 template <typename T>
234 struct FitsType;
235 
236 template <>
237 struct FitsType<bool> {
238  static int const CONSTANT = TLOGICAL;
239 };
240 template <>
241 struct FitsType<char> {
242  static int const CONSTANT = TSTRING;
243 };
244 template <>
245 struct FitsType<signed char> {
246  static int const CONSTANT = TSBYTE;
247 };
248 template <>
249 struct FitsType<unsigned char> {
250  static int const CONSTANT = TBYTE;
251 };
252 template <>
253 struct FitsType<short> {
254  static int const CONSTANT = TSHORT;
255 };
256 template <>
257 struct FitsType<unsigned short> {
258  static int const CONSTANT = TUSHORT;
259 };
260 template <>
261 struct FitsType<int> {
262  static int const CONSTANT = TINT;
263 };
264 template <>
265 struct FitsType<unsigned int> {
266  static int const CONSTANT = TUINT;
267 };
268 template <>
269 struct FitsType<long> {
270  static int const CONSTANT = TLONG;
271 };
272 template <>
273 struct FitsType<unsigned long> {
274  static int const CONSTANT = TULONG;
275 };
276 template <>
277 struct FitsType<long long> {
278  static int const CONSTANT = TLONGLONG;
279 };
280 template <>
281 struct FitsType<unsigned long long> {
282  static int const CONSTANT = TLONGLONG;
283 };
284 template <>
285 struct FitsType<float> {
286  static int const CONSTANT = TFLOAT;
287 };
288 template <>
289 struct FitsType<double> {
290  static int const CONSTANT = TDOUBLE;
291 };
292 template <>
293 struct FitsType<lsst::geom::Angle> {
294  static int const CONSTANT = TDOUBLE;
295 };
296 template <>
297 struct FitsType<std::complex<float> > {
298  static int const CONSTANT = TCOMPLEX;
299 };
300 template <>
301 struct FitsType<std::complex<double> > {
302  static int const CONSTANT = TDBLCOMPLEX;
303 };
304 
305 // We use TBIT when writing booleans to table cells, but TLOGICAL in headers.
306 template <typename T>
307 struct FitsTableType : public FitsType<T> {};
308 template <>
309 struct FitsTableType<bool> {
310  static int const CONSTANT = TBIT;
311 };
312 
313 template <typename T>
314 struct FitsBitPix;
315 
316 template <>
317 struct FitsBitPix<unsigned char> {
318  static int const CONSTANT = BYTE_IMG;
319 };
320 template <>
321 struct FitsBitPix<short> {
322  static int const CONSTANT = SHORT_IMG;
323 };
324 template <>
325 struct FitsBitPix<unsigned short> {
326  static int const CONSTANT = USHORT_IMG;
327 };
328 template <>
329 struct FitsBitPix<int> {
330  static int const CONSTANT = LONG_IMG;
331 }; // not a typo!
332 template <>
333 struct FitsBitPix<unsigned int> {
334  static int const CONSTANT = ULONG_IMG;
335 };
336 template <>
337 struct FitsBitPix<std::int64_t> {
338  static int const CONSTANT = LONGLONG_IMG;
339 };
340 template <>
341 struct FitsBitPix<std::uint64_t> {
342  static int const CONSTANT = LONGLONG_IMG;
343 };
344 template <>
345 struct FitsBitPix<float> {
346  static int const CONSTANT = FLOAT_IMG;
347 };
348 template <>
349 struct FitsBitPix<double> {
350  static int const CONSTANT = DOUBLE_IMG;
351 };
352 
353 bool isFitsImageTypeSigned(int constant) {
354  switch (constant) {
355  case BYTE_IMG: return false;
356  case SHORT_IMG: return true;
357  case USHORT_IMG: return false;
358  case LONG_IMG: return true;
359  case ULONG_IMG: return false;
360  case LONGLONG_IMG: return true;
361  case FLOAT_IMG: return true;
362  case DOUBLE_IMG: return true;
363  }
364  throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, "Invalid constant.");
365 }
366 
367 static bool allowImageCompression = true;
368 
369 int fitsTypeForBitpix(int bitpix) {
370  switch (bitpix) {
371  case 8:
372  return TBYTE;
373  case 16:
374  return TSHORT;
375  case 32:
376  return TINT;
377  case 64:
378  return TLONGLONG;
379  case -32:
380  return TFLOAT;
381  case -64:
382  return TDOUBLE;
383  default:
385  os << "Invalid bitpix value: " << bitpix;
386  throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, os.str());
387  }
388 }
389 
390 /*
391  * Information about one item of metadata: is a comment? is valid?
392  *
393  * See isCommentIsValid for more information.
394  */
395 struct ItemInfo {
396  ItemInfo(bool isComment, bool isValid) : isComment(isComment), isValid(isValid) {}
397  bool isComment;
398  bool isValid;
399 };
400 
401 /*
402  * Is an item a commnt (or history) and is it usable in a FITS header?
403  *
404  * For an item to be valid:
405  * - If name is COMMENT or HISTORY then the item must be of type std::string
406  * - All other items are always valid
407  */
408 ItemInfo isCommentIsValid(daf::base::PropertyList const &pl, std::string const &name) {
409  if (!pl.exists(name)) {
410  return ItemInfo(false, false);
411  }
412  std::type_info const &type = pl.typeOf(name);
413  if ((name == "COMMENT") || (name == "HISTORY")) {
414  return ItemInfo(true, type == typeid(std::string));
415  }
416  return ItemInfo(false, true);
417 }
418 
419 } // namespace
420 
421 // ----------------------------------------------------------------------------------------------------------
422 // ---- Implementations for stuff in fits.h -----------------------------------------------------------------
423 // ----------------------------------------------------------------------------------------------------------
424 
425 std::string makeErrorMessage(std::string const &fileName, int status, std::string const &msg) {
427  os << "cfitsio error";
428  if (fileName != "") {
429  os << " (" << fileName << ")";
430  }
431  if (status != 0) {
432  char fitsErrMsg[FLEN_ERRMSG];
433  fits_get_errstatus(status, fitsErrMsg);
434  os << ": " << fitsErrMsg << " (" << status << ")";
435  }
436  if (msg != "") {
437  os << " : " << msg;
438  }
439  os << "\ncfitsio error stack:\n";
440  char cfitsioMsg[FLEN_ERRMSG];
441  while (fits_read_errmsg(cfitsioMsg) != 0) {
442  os << " " << cfitsioMsg << "\n";
443  }
444  return os.str();
445 }
446 
447 std::string makeErrorMessage(void *fptr, int status, std::string const &msg) {
448  std::string fileName = "";
449  fitsfile *fd = reinterpret_cast<fitsfile *>(fptr);
450  if (fd != 0 && fd->Fptr != 0 && fd->Fptr->filename != 0) {
451  fileName = fd->Fptr->filename;
452  }
453  return makeErrorMessage(fileName, status, msg);
454 }
455 
457  std::set<std::string> const &excludeNames) {
458  daf::base::PropertyList const *pl = dynamic_cast<daf::base::PropertyList const *>(&metadata);
459  std::vector<std::string> allParamNames;
460  if (pl) {
461  allParamNames = pl->getOrderedNames();
462  } else {
463  allParamNames = metadata.paramNames(false);
464  }
465  std::vector<std::string> desiredParamNames;
466  for (auto const &name : allParamNames) {
467  if (excludeNames.count(name) == 0) {
468  desiredParamNames.push_back(name);
469  }
470  }
471  return makeLimitedFitsHeaderImpl(desiredParamNames, metadata);
472 }
473 
475  if (_managed) std::free(_ptr);
476  _ptr = 0;
477  _len = 0;
478  _managed = true;
479 }
480 
482  reset();
483  _ptr = std::malloc(len);
484  _len = len;
485  _managed = true;
486 }
487 
488 template <typename T>
489 int getBitPix() {
490  return FitsBitPix<T>::CONSTANT;
491 }
492 
493 // ----------------------------------------------------------------------------------------------------------
494 // ---- Implementations for Fits class ----------------------------------------------------------------------
495 // ----------------------------------------------------------------------------------------------------------
496 
498  std::string fileName = "<unknown>";
499  fitsfile *fd = reinterpret_cast<fitsfile *>(fptr);
500  if (fd != 0 && fd->Fptr != 0 && fd->Fptr->filename != 0) {
501  fileName = fd->Fptr->filename;
502  }
503  return fileName;
504 }
505 
507  int n = 1;
508  fits_get_hdu_num(reinterpret_cast<fitsfile *>(fptr), &n);
509  return n - 1;
510 }
511 
512 void Fits::setHdu(int hdu, bool relative) {
513  if (relative) {
514  fits_movrel_hdu(reinterpret_cast<fitsfile *>(fptr), hdu, 0, &status);
515  if (behavior & AUTO_CHECK) {
516  LSST_FITS_CHECK_STATUS(*this, boost::format("Incrementing HDU by %d") % hdu);
517  }
518  } else {
519  if (hdu != DEFAULT_HDU) {
520  fits_movabs_hdu(reinterpret_cast<fitsfile *>(fptr), hdu + 1, 0, &status);
521  }
522  if (hdu == DEFAULT_HDU && getHdu() == 0 && getImageDim() == 0) {
523  // want a silent failure here
524  int tmpStatus = status;
525  fits_movrel_hdu(reinterpret_cast<fitsfile *>(fptr), 1, 0, &tmpStatus);
526  }
527  if (behavior & AUTO_CHECK) {
528  LSST_FITS_CHECK_STATUS(*this, boost::format("Moving to HDU %d") % hdu);
529  }
530  }
531 }
532 
534  int n = 0;
535  fits_get_num_hdus(reinterpret_cast<fitsfile *>(fptr), &n, &status);
536  if (behavior & AUTO_CHECK) {
537  LSST_FITS_CHECK_STATUS(*this, "Getting number of HDUs in file.");
538  }
539  return n;
540 }
541 
542 // ---- Writing and updating header keys --------------------------------------------------------------------
543 
544 namespace {
545 
546 // Impl functions in the anonymous namespace do special handling for strings, bools, and IEEE fp values.
547 
554 std::string nonFiniteDoubleToString(double value) {
555  if (std::isfinite(value)) {
556  return "";
557  }
558  if (std::isnan(value)) {
559  return "NAN";
560  }
561  if (value < 0) {
562  return "-INFINITY";
563  }
564  return "+INFINITY";
565 }
566 
572 double stringToNonFiniteDouble(std::string const &value) {
573  if (value == "NAN") {
575  }
576  if (value == "+INFINITY") {
578  }
579  if (value == "-INFINITY") {
581  }
582  return 0;
583 }
584 
585 template <typename T>
586 void updateKeyImpl(Fits &fits, char const *key, T const &value, char const *comment) {
587  fits_update_key(reinterpret_cast<fitsfile *>(fits.fptr), FitsType<T>::CONSTANT, const_cast<char *>(key),
588  const_cast<T *>(&value), const_cast<char *>(comment), &fits.status);
589 }
590 
591 void updateKeyImpl(Fits &fits, char const *key, std::string const &value, char const *comment) {
592  fits_update_key_longstr(reinterpret_cast<fitsfile *>(fits.fptr), const_cast<char *>(key),
593  const_cast<char *>(value.c_str()), const_cast<char *>(comment), &fits.status);
594 }
595 
596 void updateKeyImpl(Fits &fits, char const *key, bool const &value, char const *comment) {
597  int v = value;
598  fits_update_key(reinterpret_cast<fitsfile *>(fits.fptr), TLOGICAL, const_cast<char *>(key), &v,
599  const_cast<char *>(comment), &fits.status);
600 }
601 
602 void updateKeyImpl(Fits &fits, char const *key, double const &value, char const *comment) {
603  std::string strValue = nonFiniteDoubleToString(value);
604  if (!strValue.empty()) {
605  updateKeyImpl(fits, key, strValue, comment);
606  } else {
607  fits_update_key(reinterpret_cast<fitsfile *>(fits.fptr), FitsType<double>::CONSTANT,
608  const_cast<char *>(key), const_cast<double *>(&value), const_cast<char *>(comment),
609  &fits.status);
610  }
611 }
612 
613 template <typename T>
614 void writeKeyImpl(Fits &fits, char const *key, T const &value, char const *comment) {
615  fits_write_key(reinterpret_cast<fitsfile *>(fits.fptr), FitsType<T>::CONSTANT, const_cast<char *>(key),
616  const_cast<T *>(&value), const_cast<char *>(comment), &fits.status);
617 }
618 
619 void writeKeyImpl(Fits &fits, char const *key, char const *comment) {
620  // Write a key with an undefined value
621  fits_write_key_null(reinterpret_cast<fitsfile *>(fits.fptr), const_cast<char *>(key),
622  const_cast<char *>(comment), &fits.status);
623 }
624 
625 void writeKeyImpl(Fits &fits, char const *key, std::string const &value, char const *comment) {
626  if (strncmp(key, "COMMENT", 7) == 0) {
627  fits_write_comment(reinterpret_cast<fitsfile *>(fits.fptr), const_cast<char *>(value.c_str()),
628  &fits.status);
629  } else if (strncmp(key, "HISTORY", 7) == 0) {
630  fits_write_history(reinterpret_cast<fitsfile *>(fits.fptr), const_cast<char *>(value.c_str()),
631  &fits.status);
632  } else {
633  fits_write_key_longstr(reinterpret_cast<fitsfile *>(fits.fptr), const_cast<char *>(key),
634  const_cast<char *>(value.c_str()), const_cast<char *>(comment), &fits.status);
635  }
636 }
637 
638 void writeKeyImpl(Fits &fits, char const *key, bool const &value, char const *comment) {
639  int v = value;
640  fits_write_key(reinterpret_cast<fitsfile *>(fits.fptr), TLOGICAL, const_cast<char *>(key), &v,
641  const_cast<char *>(comment), &fits.status);
642 }
643 
644 void writeKeyImpl(Fits &fits, char const *key, double const &value, char const *comment) {
645  std::string strValue = nonFiniteDoubleToString(value);
646  if (!strValue.empty()) {
647  writeKeyImpl(fits, key, strValue, comment);
648  } else {
649  fits_write_key(reinterpret_cast<fitsfile *>(fits.fptr), FitsType<double>::CONSTANT,
650  const_cast<char *>(key), const_cast<double *>(&value), const_cast<char *>(comment),
651  &fits.status);
652  }
653 }
654 
655 } // namespace
656 
657 template <typename T>
658 void Fits::updateKey(std::string const &key, T const &value, std::string const &comment) {
659  updateKeyImpl(*this, key.c_str(), value, comment.c_str());
660  if (behavior & AUTO_CHECK) {
661  LSST_FITS_CHECK_STATUS(*this, boost::format("Updating key '%s': '%s'") % key % value);
662  }
663 }
664 
665 template <typename T>
666 void Fits::writeKey(std::string const &key, T const &value, std::string const &comment) {
667  writeKeyImpl(*this, key.c_str(), value, comment.c_str());
668  if (behavior & AUTO_CHECK) {
669  LSST_FITS_CHECK_STATUS(*this, boost::format("Writing key '%s': '%s'") % key % value);
670  }
671 }
672 
673 template <typename T>
674 void Fits::updateKey(std::string const &key, T const &value) {
675  updateKeyImpl(*this, key.c_str(), value, 0);
676  if (behavior & AUTO_CHECK) {
677  LSST_FITS_CHECK_STATUS(*this, boost::format("Updating key '%s': '%s'") % key % value);
678  }
679 }
680 
681 template <typename T>
682 void Fits::writeKey(std::string const &key, T const &value) {
683  writeKeyImpl(*this, key.c_str(), value, 0);
684  if (behavior & AUTO_CHECK) {
685  LSST_FITS_CHECK_STATUS(*this, boost::format("Writing key '%s': '%s'") % key % value);
686  }
687 }
688 
689 template <typename T>
690 void Fits::updateColumnKey(std::string const &prefix, int n, T const &value, std::string const &comment) {
691  updateKey((boost::format("%s%d") % prefix % (n + 1)).str(), value, comment);
692  if (behavior & AUTO_CHECK) {
693  LSST_FITS_CHECK_STATUS(*this, boost::format("Updating key '%s%d': '%s'") % prefix % (n + 1) % value);
694  }
695 }
696 
697 template <typename T>
698 void Fits::writeColumnKey(std::string const &prefix, int n, T const &value, std::string const &comment) {
699  writeKey((boost::format("%s%d") % prefix % (n + 1)).str(), value, comment);
700  if (behavior & AUTO_CHECK) {
701  LSST_FITS_CHECK_STATUS(*this, boost::format("Writing key '%s%d': '%s'") % prefix % (n + 1) % value);
702  }
703 }
704 
705 template <typename T>
706 void Fits::updateColumnKey(std::string const &prefix, int n, T const &value) {
707  updateKey((boost::format("%s%d") % prefix % (n + 1)).str(), value);
708  if (behavior & AUTO_CHECK) {
709  LSST_FITS_CHECK_STATUS(*this, boost::format("Updating key '%s%d': '%s'") % prefix % (n + 1) % value);
710  }
711 }
712 
713 template <typename T>
714 void Fits::writeColumnKey(std::string const &prefix, int n, T const &value) {
715  writeKey((boost::format("%s%d") % prefix % (n + 1)).str(), value);
716  if (behavior & AUTO_CHECK) {
717  LSST_FITS_CHECK_STATUS(*this, boost::format("Writing key '%s%d': '%s'") % prefix % (n + 1) % value);
718  }
719 }
720 
721 // ---- Reading header keys ---------------------------------------------------------------------------------
722 
723 namespace {
724 
725 template <typename T>
726 void readKeyImpl(Fits &fits, char const *key, T &value) {
727  fits_read_key(reinterpret_cast<fitsfile *>(fits.fptr), FitsType<T>::CONSTANT, const_cast<char *>(key),
728  &value, 0, &fits.status);
729 }
730 
731 void readKeyImpl(Fits &fits, char const *key, std::string &value) {
732  char *buf = 0;
733  fits_read_key_longstr(reinterpret_cast<fitsfile *>(fits.fptr), const_cast<char *>(key), &buf, 0,
734  &fits.status);
735  if (buf) {
736  value = strip(buf);
737  free(buf);
738  }
739 }
740 
741 void readKeyImpl(Fits &fits, char const *key, double &value) {
742  // We need to check for the possibility that the value is a special string (for NAN, +/-Inf).
743  // If a quote mark (') is present then it's a string.
744 
745  char buf[FLEN_VALUE];
746  fits_read_keyword(reinterpret_cast<fitsfile *>(fits.fptr), const_cast<char *>(key), buf, 0, &fits.status);
747  if (fits.status != 0) {
748  return;
749  }
750  if (std::string(buf).find('\'') != std::string::npos) {
751  std::string unquoted;
752  readKeyImpl(fits, key, unquoted); // Let someone else remove quotes and whitespace
753  if (fits.status != 0) {
754  return;
755  }
756  value = stringToNonFiniteDouble(unquoted);
757  if (value == 0) {
758  throw LSST_EXCEPT(
760  (boost::format("Unrecognised string value for keyword '%s' when parsing as double: %s") %
761  key % unquoted)
762  .str());
763  }
764  } else {
765  fits_read_key(reinterpret_cast<fitsfile *>(fits.fptr), FitsType<double>::CONSTANT,
766  const_cast<char *>(key), &value, 0, &fits.status);
767  }
768 }
769 
770 } // namespace
771 
772 template <typename T>
773 void Fits::readKey(std::string const &key, T &value) {
774  readKeyImpl(*this, key.c_str(), value);
775  if (behavior & AUTO_CHECK) {
776  LSST_FITS_CHECK_STATUS(*this, boost::format("Reading key '%s'") % key);
777  }
778 }
779 
781  char key[81]; // allow for terminating NUL
782  char value[81];
783  char comment[81];
784  int nKeys = 0;
785  fits_get_hdrspace(reinterpret_cast<fitsfile *>(fptr), &nKeys, 0, &status);
786  std::string keyStr;
787  std::string valueStr;
788  std::string commentStr;
789  int i = 1;
790  while (i <= nKeys) {
791  fits_read_keyn(reinterpret_cast<fitsfile *>(fptr), i, key, value, comment, &status);
792  keyStr = key;
793  valueStr = value;
794  commentStr = comment;
795  ++i;
796  while (valueStr.size() > 2 && valueStr[valueStr.size() - 2] == '&' && i <= nKeys) {
797  // we're using key to hold the entire record here; the actual key is safe in keyStr
798  fits_read_record(reinterpret_cast<fitsfile *>(fptr), i, key, &status);
799  if (strncmp(key, "CONTINUE", 8) != 0) {
800  // require both trailing '&' and CONTINUE to invoke long-string handling
801  break;
802  }
803  std::string card = key;
804  valueStr.erase(valueStr.size() - 2);
805  std::size_t firstQuote = card.find('\'');
806  if (firstQuote == std::string::npos) {
807  throw LSST_EXCEPT(
808  FitsError,
810  fptr, status,
811  boost::format("Invalid CONTINUE at header key %d: \"%s\".") % i % card));
812  }
813  std::size_t lastQuote = card.find('\'', firstQuote + 1);
814  if (lastQuote == std::string::npos) {
815  throw LSST_EXCEPT(
816  FitsError,
818  fptr, status,
819  boost::format("Invalid CONTINUE at header key %d: \"%s\".") % i % card));
820  }
821  valueStr += card.substr(firstQuote + 1, lastQuote - firstQuote);
822  std::size_t slash = card.find('/', lastQuote + 1);
823  if (slash != std::string::npos) {
824  commentStr += strip(card.substr(slash + 1));
825  }
826  ++i;
827  }
828  if (behavior & AUTO_CHECK) {
829  LSST_FITS_CHECK_STATUS(*this, boost::format("Reading key '%s'") % keyStr);
830  }
831  functor(keyStr, valueStr, commentStr);
832  }
833 }
834 
835 // ---- Reading and writing PropertySet/PropertyList --------------------------------------------------------
836 
837 namespace {
838 
839 bool isKeyIgnored(std::string const &key, bool write = false) {
840  return ((ignoreKeys.find(key) != ignoreKeys.end()) || ignoreKeyStarts.matches(key) ||
841  (write && ignoreKeyStartsWrite.matches(key)));
842 }
843 
844 class MetadataIterationFunctor : public HeaderIterationFunctor {
845 public:
846  void operator()(std::string const &key, std::string const &value, std::string const &comment) override;
847 
848  template <typename T>
849  void add(std::string const &key, T value, std::string const &comment) {
850  // PropertyList/Set can not support array items where some elements are
851  // defined and some undefined. If we are adding defined value where
852  // previously we have an undefined value we must use set instead.
853  if (list) {
854  if (list->exists(key) && list->isUndefined(key)) {
855  LOGLS_WARN("afw.fits",
856  boost::format("In %s, replacing undefined value for key '%s'.") %
857  BOOST_CURRENT_FUNCTION % key);
858  list->set(key, value, comment);
859  } else {
860  list->add(key, value, comment);
861  }
862  } else {
863  if (set->exists(key) && set->isUndefined(key)) {
864  LOGLS_WARN("afw.fits",
865  boost::format("In %s, replacing undefined value for key '%s'.") %
866  BOOST_CURRENT_FUNCTION % key);
867  set->set(key, value);
868  } else {
869  set->add(key, value);
870  }
871  }
872  }
873 
874  void add(std::string const &key, std::string const &comment) {
875  // If this undefined value is adding to a pre-existing key that has
876  // a defined value we must skip the add so as not to break
877  // PropertyList/Set.
878  if (list) {
879  if (list->exists(key) && !list->isUndefined(key)) {
880  // Do nothing. Assume the previously defined value takes
881  // precedence.
882  LOGLS_WARN("afw.fits",
883  boost::format("In %s, dropping undefined value for key '%s'.") %
884  BOOST_CURRENT_FUNCTION % key);
885  } else {
886  list->add(key, nullptr, comment);
887  }
888  } else {
889  if (set->exists(key) && !set->isUndefined(key)) {
890  // Do nothing. Assume the previously defined value takes
891  // precedence.
892  LOGLS_WARN("afw.fits",
893  boost::format("In %s, dropping undefined value for key '%s'.") %
894  BOOST_CURRENT_FUNCTION % key);
895  } else {
896  set->add(key, nullptr);
897  }
898  }
899  }
900 
901  bool strip;
902  daf::base::PropertySet *set;
903  daf::base::PropertyList *list;
904 };
905 
906 void MetadataIterationFunctor::operator()(std::string const &key, std::string const &value,
907  std::string const &comment) {
908  static boost::regex const boolRegex("[tTfF]");
909  static boost::regex const intRegex("[+-]?[0-9]+");
910  static boost::regex const doubleRegex("[+-]?([0-9]*\\.?[0-9]+|[0-9]+\\.?[0-9]*)([eE][+-]?[0-9]+)?");
911  static boost::regex const fitsStringRegex("'(.*?) *'");
912  // regex for two-line comment added to all FITS headers by CFITSIO
913  static boost::regex const fitsDefinitionCommentRegex(
914  " *(FITS \\(Flexible Image Transport System\\)|and Astrophysics', volume 376, page 359).*");
915  boost::smatch matchStrings;
916 
917  if (strip && isKeyIgnored(key)) {
918  return;
919  }
920 
921  std::istringstream converter(value);
922  if (boost::regex_match(value, boolRegex)) {
923  // convert the string to an bool
924  add(key, bool(value == "T" || value == "t"), comment);
925  } else if (boost::regex_match(value, intRegex)) {
926  // convert the string to an int
928  converter >> val;
929  if (val < (1LL << 31) && val > -(1LL << 31)) {
930  add(key, static_cast<int>(val), comment);
931  } else {
932  add(key, val, comment);
933  }
934  } else if (boost::regex_match(value, doubleRegex)) {
935  // convert the string to a double
936  double val;
937  converter >> val;
938  add(key, val, comment);
939  } else if (boost::regex_match(value, matchStrings, fitsStringRegex)) {
940  std::string const str = matchStrings[1].str(); // strip off the enclosing single quotes
941  double val = stringToNonFiniteDouble(str);
942  if (val != 0.0) {
943  add(key, val, comment);
944  } else {
945  add(key, str, comment);
946  }
947  } else if (key == "HISTORY") {
948  add(key, comment, "");
949  } else if (key == "COMMENT" && !(strip && boost::regex_match(comment, fitsDefinitionCommentRegex))) {
950  add(key, comment, "");
951  } else if (key.empty() && value.empty()) {
952  // This is a blank keyword comment. Since comments do not retain
953  // their position on read there is nothing to be gained by storing
954  // this in the PropertyList as a blank keyword. Therefore store
955  // them with the other comments.
956  add("COMMENT", comment, "");
957  } else if (value.empty()) {
958  // do nothing for empty values that are comments
959  // Otherwise write null value to PropertySet
960  if (key != "COMMENT") {
961  add(key, comment);
962  }
963  } else {
964  throw LSST_EXCEPT(
965  afw::fits::FitsError,
966  (boost::format("Could not parse header value for key '%s': '%s'") % key % value).str());
967  }
968 }
969 
970 void writeKeyFromProperty(Fits &fits, daf::base::PropertySet const &metadata, std::string const &key,
971  char const *comment = 0) {
972  std::type_info const &valueType = metadata.typeOf(key);
973  if (valueType == typeid(bool)) {
974  if (metadata.isArray(key)) {
975  std::vector<bool> tmp = metadata.getArray<bool>(key);
976  // work around unfortunate specialness of std::vector<bool>
977  for (std::size_t i = 0; i != tmp.size(); ++i) {
978  writeKeyImpl(fits, key.c_str(), static_cast<bool>(tmp[i]), comment);
979  }
980  } else {
981  writeKeyImpl(fits, key.c_str(), metadata.get<bool>(key), comment);
982  }
983  } else if (valueType == typeid(std::uint8_t)) {
984  if (metadata.isArray(key)) {
985  std::vector<std::uint8_t> tmp = metadata.getArray<std::uint8_t>(key);
986  for (std::size_t i = 0; i != tmp.size(); ++i) {
987  writeKeyImpl(fits, key.c_str(), tmp[i], comment);
988  }
989  } else {
990  writeKeyImpl(fits, key.c_str(), metadata.get<std::uint8_t>(key), comment);
991  }
992  } else if (valueType == typeid(int)) {
993  if (metadata.isArray(key)) {
994  std::vector<int> tmp = metadata.getArray<int>(key);
995  for (std::size_t i = 0; i != tmp.size(); ++i) {
996  writeKeyImpl(fits, key.c_str(), tmp[i], comment);
997  }
998  } else {
999  writeKeyImpl(fits, key.c_str(), metadata.get<int>(key), comment);
1000  }
1001  } else if (valueType == typeid(long)) {
1002  if (metadata.isArray(key)) {
1003  std::vector<long> tmp = metadata.getArray<long>(key);
1004  for (std::size_t i = 0; i != tmp.size(); ++i) {
1005  writeKeyImpl(fits, key.c_str(), tmp[i], comment);
1006  }
1007  } else {
1008  writeKeyImpl(fits, key.c_str(), metadata.get<long>(key), comment);
1009  }
1010  } else if (valueType == typeid(long long)) {
1011  if (metadata.isArray(key)) {
1012  std::vector<long long> tmp = metadata.getArray<long long>(key);
1013  for (std::size_t i = 0; i != tmp.size(); ++i) {
1014  writeKeyImpl(fits, key.c_str(), tmp[i], comment);
1015  }
1016  } else {
1017  writeKeyImpl(fits, key.c_str(), metadata.get<long long>(key), comment);
1018  }
1019  } else if (valueType == typeid(std::int64_t)) {
1020  if (metadata.isArray(key)) {
1021  std::vector<std::int64_t> tmp = metadata.getArray<std::int64_t>(key);
1022  for (std::size_t i = 0; i != tmp.size(); ++i) {
1023  writeKeyImpl(fits, key.c_str(), tmp[i], comment);
1024  }
1025  } else {
1026  writeKeyImpl(fits, key.c_str(), metadata.get<std::int64_t>(key), comment);
1027  }
1028  } else if (valueType == typeid(double)) {
1029  if (metadata.isArray(key)) {
1030  std::vector<double> tmp = metadata.getArray<double>(key);
1031  for (std::size_t i = 0; i != tmp.size(); ++i) {
1032  writeKeyImpl(fits, key.c_str(), tmp[i], comment);
1033  }
1034  } else {
1035  writeKeyImpl(fits, key.c_str(), metadata.get<double>(key), comment);
1036  }
1037  } else if (valueType == typeid(std::string)) {
1038  if (metadata.isArray(key)) {
1039  std::vector<std::string> tmp = metadata.getArray<std::string>(key);
1040  for (std::size_t i = 0; i != tmp.size(); ++i) {
1041  writeKeyImpl(fits, key.c_str(), tmp[i], comment);
1042  }
1043  } else {
1044  writeKeyImpl(fits, key.c_str(), metadata.get<std::string>(key), comment);
1045  }
1046  } else if (valueType == typeid(nullptr_t)) {
1047  if (metadata.isArray(key)) {
1048  // Write multiple undefined values for the same key
1049  std::vector<std::string> tmp = metadata.getArray<std::string>(key);
1050  for (std::size_t i = 0; i != tmp.size(); ++i) {
1051  writeKeyImpl(fits, key.c_str(), comment);
1052  }
1053  } else {
1054  writeKeyImpl(fits, key.c_str(), comment);
1055  }
1056  } else {
1057  // FIXME: inherited this error handling from fitsIo.cc; need a better option.
1058  LOGLS_WARN("afw.writeKeyFromProperty",
1059  makeErrorMessage(fits.fptr, fits.status,
1060  boost::format("In %s, unknown type '%s' for key '%s'.") %
1061  BOOST_CURRENT_FUNCTION % valueType.name() % key));
1062  }
1063  if (fits.behavior & Fits::AUTO_CHECK) {
1064  LSST_FITS_CHECK_STATUS(fits, boost::format("Writing key '%s'") % key);
1065  }
1066 }
1067 
1068 } // namespace
1069 
1071  MetadataIterationFunctor f;
1072  f.strip = strip;
1073  f.set = &metadata;
1074  f.list = dynamic_cast<daf::base::PropertyList *>(&metadata);
1075  forEachKey(f);
1076 }
1077 
1079  typedef std::vector<std::string> NameList;
1080  daf::base::PropertyList const *pl = dynamic_cast<daf::base::PropertyList const *>(&metadata);
1081  NameList paramNames;
1082  if (pl) {
1083  paramNames = pl->getOrderedNames();
1084  } else {
1085  paramNames = metadata.paramNames(false);
1086  }
1087  for (NameList::const_iterator i = paramNames.begin(); i != paramNames.end(); ++i) {
1088  if (!isKeyIgnored(*i, true)) {
1089  if (pl) {
1090  writeKeyFromProperty(*this, metadata, *i, pl->getComment(*i).c_str());
1091  } else {
1092  writeKeyFromProperty(*this, metadata, *i);
1093  }
1094  }
1095  }
1096 }
1097 
1098 // ---- Manipulating tables ---------------------------------------------------------------------------------
1099 
1101  char *ttype = 0;
1102  char *tform = 0;
1103  fits_create_tbl(reinterpret_cast<fitsfile *>(fptr), BINARY_TBL, 0, 0, &ttype, &tform, 0, 0, &status);
1104  if (behavior & AUTO_CHECK) {
1105  LSST_FITS_CHECK_STATUS(*this, "Creating binary table");
1106  }
1107 }
1108 
1109 template <typename T>
1110 int Fits::addColumn(std::string const &ttype, int size) {
1111  int nCols = 0;
1112  fits_get_num_cols(reinterpret_cast<fitsfile *>(fptr), &nCols, &status);
1113  std::string tform = makeColumnFormat<T>(size);
1114  fits_insert_col(reinterpret_cast<fitsfile *>(fptr), nCols + 1, const_cast<char *>(ttype.c_str()),
1115  const_cast<char *>(tform.c_str()), &status);
1116  if (behavior & AUTO_CHECK) {
1117  LSST_FITS_CHECK_STATUS(*this, boost::format("Adding column '%s' with size %d") % ttype % size);
1118  }
1119  return nCols;
1120 }
1121 
1122 template <typename T>
1123 int Fits::addColumn(std::string const &ttype, int size, std::string const &comment) {
1124  int nCols = addColumn<T>(ttype, size);
1125  updateColumnKey("TTYPE", nCols, ttype, comment);
1126  if (behavior & AUTO_CHECK) {
1127  LSST_FITS_CHECK_STATUS(*this, boost::format("Adding column '%s' with size %d") % ttype % size);
1128  }
1129  return nCols;
1130 }
1131 
1133  long first = 0;
1134  fits_get_num_rows(reinterpret_cast<fitsfile *>(fptr), &first, &status);
1135  fits_insert_rows(reinterpret_cast<fitsfile *>(fptr), first, nRows, &status);
1136  if (behavior & AUTO_CHECK) {
1137  LSST_FITS_CHECK_STATUS(*this, boost::format("Adding %d rows to binary table") % nRows);
1138  }
1139  return first;
1140 }
1141 
1143  long r = 0;
1144  fits_get_num_rows(reinterpret_cast<fitsfile *>(fptr), &r, &status);
1145  if (behavior & AUTO_CHECK) {
1146  LSST_FITS_CHECK_STATUS(*this, "Checking how many rows are in table");
1147  }
1148  return r;
1149 }
1150 
1151 template <typename T>
1152 void Fits::writeTableArray(std::size_t row, int col, int nElements, T const *value) {
1153  fits_write_col(reinterpret_cast<fitsfile *>(fptr), FitsTableType<T>::CONSTANT, col + 1, row + 1, 1,
1154  nElements, const_cast<T *>(value), &status);
1155  if (behavior & AUTO_CHECK) {
1156  LSST_FITS_CHECK_STATUS(*this, boost::format("Writing %d-element array at table cell (%d, %d)") %
1157  nElements % row % col);
1158  }
1159 }
1160 
1162  // cfitsio doesn't let us specify the size of a string, it just looks for null terminator.
1163  // Using std::string::c_str() guarantees that we have one. But we can't store arbitrary
1164  // data in a string field because cfitsio will also chop off anything after the first null
1165  // terminator.
1166  char const *tmp = value.c_str();
1167  fits_write_col(reinterpret_cast<fitsfile *>(fptr), TSTRING, col + 1, row + 1, 1, 1,
1168  const_cast<char const **>(&tmp), &status);
1169  if (behavior & AUTO_CHECK) {
1170  LSST_FITS_CHECK_STATUS(*this, boost::format("Writing value at table cell (%d, %d)") % row % col);
1171  }
1172 }
1173 
1174 template <typename T>
1175 void Fits::readTableArray(std::size_t row, int col, int nElements, T *value) {
1176  int anynul = false;
1177  fits_read_col(reinterpret_cast<fitsfile *>(fptr), FitsTableType<T>::CONSTANT, col + 1, row + 1, 1,
1178  nElements, 0, value, &anynul, &status);
1179  if (behavior & AUTO_CHECK) {
1180  LSST_FITS_CHECK_STATUS(*this, boost::format("Reading value at table cell (%d, %d)") % row % col);
1181  }
1182 }
1183 
1184 void Fits::readTableScalar(std::size_t row, int col, std::string &value, bool isVariableLength) {
1185  int anynul = false;
1186  long size = isVariableLength ? getTableArraySize(row, col) : getTableArraySize(col);
1187  // We can't directly write into a std::string until C++17.
1188  std::vector<char> buf(size + 1, 0);
1189  // cfitsio wants a char** because they imagine we might want an array of strings,
1190  // but we only want one element.
1191  char *tmp = &buf.front();
1192  fits_read_col(reinterpret_cast<fitsfile *>(fptr), TSTRING, col + 1, row + 1, 1, 1, 0, &tmp, &anynul,
1193  &status);
1194  if (behavior & AUTO_CHECK) {
1195  LSST_FITS_CHECK_STATUS(*this, boost::format("Reading value at table cell (%d, %d)") % row % col);
1196  }
1197  value = std::string(tmp);
1198 }
1199 
1201  int typecode = 0;
1202  long result = 0;
1203  long width = 0;
1204  fits_get_coltype(reinterpret_cast<fitsfile *>(fptr), col + 1, &typecode, &result, &width, &status);
1205  if (behavior & AUTO_CHECK) {
1206  LSST_FITS_CHECK_STATUS(*this, boost::format("Looking up array size for column %d") % col);
1207  }
1208  return result;
1209 }
1210 
1212  long result = 0;
1213  long offset = 0;
1214  fits_read_descript(reinterpret_cast<fitsfile *>(fptr), col + 1, row + 1, &result, &offset, &status);
1215  if (behavior & AUTO_CHECK) {
1216  LSST_FITS_CHECK_STATUS(*this, boost::format("Looking up array size for cell (%d, %d)") % row % col);
1217  }
1218  return result;
1219 }
1220 
1221 // ---- Manipulating images ---------------------------------------------------------------------------------
1222 
1224  long naxes = 0;
1225  fits_create_img(reinterpret_cast<fitsfile *>(fptr), 8, 0, &naxes, &status);
1226  if (behavior & AUTO_CHECK) {
1227  LSST_FITS_CHECK_STATUS(*this, "Creating empty image HDU");
1228  }
1229 }
1230 
1231 void Fits::createImageImpl(int bitpix, int naxis, long const *naxes) {
1232  fits_create_img(reinterpret_cast<fitsfile *>(fptr), bitpix, naxis, const_cast<long *>(naxes), &status);
1233  if (behavior & AUTO_CHECK) {
1234  LSST_FITS_CHECK_STATUS(*this, "Creating new image HDU");
1235  }
1236 }
1237 
1238 template <typename T>
1239 void Fits::writeImageImpl(T const *data, int nElements) {
1240  fits_write_img(reinterpret_cast<fitsfile *>(fptr), FitsType<T>::CONSTANT, 1, nElements,
1241  const_cast<T *>(data), &status);
1242  if (behavior & AUTO_CHECK) {
1243  LSST_FITS_CHECK_STATUS(*this, "Writing image");
1244  }
1245 }
1246 
1247 namespace {
1248 
1253 struct ImageCompressionContext {
1254 public:
1255  ImageCompressionContext(Fits &fits_, ImageCompressionOptions const &useThis)
1256  : fits(fits_), old(fits.getImageCompression()) {
1257  fits.setImageCompression(useThis);
1258  }
1259  ~ImageCompressionContext() {
1260  int status = 0;
1261  std::swap(status, fits.status);
1262  try {
1263  fits.setImageCompression(old);
1264  } catch (...) {
1265  LOGLS_WARN("afw.fits",
1266  makeErrorMessage(fits.fptr, fits.status, "Failed to restore compression settings"));
1267  }
1268  std::swap(status, fits.status);
1269  }
1270 
1271 private:
1272  Fits &fits; // FITS file we're working on
1273  ImageCompressionOptions old; // Former compression options, to be restored
1274 };
1275 
1276 } // anonymous namespace
1277 
1278 template <typename T>
1282  auto fits = reinterpret_cast<fitsfile *>(fptr);
1283  ImageCompressionOptions const &compression =
1284  image.getBBox().getArea() > 0
1285  ? options.compression
1287  ImageCompressionOptions::NONE); // cfitsio can't compress empty images
1288  ImageCompressionContext comp(*this, compression); // RAII
1289  if (behavior & AUTO_CHECK) {
1290  LSST_FITS_CHECK_STATUS(*this, "Activating compression for write image");
1291  }
1292 
1293  ImageScale scale = options.scaling.determine(image, mask);
1294 
1295  // We need a place to put the image+header, and CFITSIO needs to know the dimenions.
1296  ndarray::Vector<long, 2> dims(image.getArray().getShape().reverse());
1297  createImageImpl(scale.bitpix == 0 ? detail::Bitpix<T>::value : scale.bitpix, 2, dims.elems);
1298 
1299  // Write the header
1302  if (header) {
1303  std::shared_ptr<daf::base::PropertySet> copy = header->deepCopy();
1304  copy->combine(wcsMetadata);
1305  header = copy;
1306  } else {
1307  header = wcsMetadata;
1308  }
1309  writeMetadata(*header);
1310 
1311  // Scale the image how we want it on disk
1312  ndarray::Array<T const, 2, 2> array = makeContiguousArray(image.getArray());
1313  auto pixels = scale.toFits(array, compression.quantizeLevel != 0, options.scaling.fuzz,
1314  options.compression.tiles, options.scaling.seed);
1315 
1316  // We only want cfitsio to do the scale and zero for unsigned 64-bit integer types. For those,
1317  // "double bzero" has sufficient precision to represent the appropriate value. We'll let
1318  // cfitsio handle it itself.
1319  // In all other cases, we will convert the image to use the appropriate scale and zero
1320  // (because we want to fuzz the numbers in the quantisation), so we don't want cfitsio
1321  // rescaling.
1323  fits_set_bscale(fits, 1.0, 0.0, &status);
1324  if (behavior & AUTO_CHECK) {
1325  LSST_FITS_CHECK_STATUS(*this, "Setting bscale,bzero");
1326  }
1327  }
1328 
1329  // Write the pixels
1330  int const fitsType = scale.bitpix == 0 ? FitsType<T>::CONSTANT : fitsTypeForBitpix(scale.bitpix);
1331  fits_write_img(fits, fitsType, 1, pixels->getNumElements(), const_cast<void *>(pixels->getData()),
1332  &status);
1333  if (behavior & AUTO_CHECK) {
1334  LSST_FITS_CHECK_STATUS(*this, "Writing image");
1335  }
1336 
1337  // Now write the headers we didn't want cfitsio to know about when we were writing the pixels
1338  // (because we don't want it using them to modify the pixels, and we don't want it overwriting
1339  // these values).
1341  std::isfinite(scale.bzero) && std::isfinite(scale.bscale) && (scale.bscale != 0.0)) {
1343  if (scale.bzero != 0.0) {
1344  fits_write_key_lng(fits, "BZERO", static_cast<long>(scale.bzero),
1345  "Scaling: MEMORY = BZERO + BSCALE * DISK", &status);
1346  }
1347  if (scale.bscale != 1.0) {
1348  fits_write_key_lng(fits, "BSCALE", static_cast<long>(scale.bscale),
1349  "Scaling: MEMORY = BZERO + BSCALE * DISK", &status);
1350  }
1351  } else {
1352  fits_write_key_dbl(fits, "BZERO", scale.bzero, 12, "Scaling: MEMORY = BZERO + BSCALE * DISK",
1353  &status);
1354  fits_write_key_dbl(fits, "BSCALE", scale.bscale, 12, "Scaling: MEMORY = BZERO + BSCALE * DISK",
1355  &status);
1356  }
1357  if (behavior & AUTO_CHECK) {
1358  LSST_FITS_CHECK_STATUS(*this, "Writing BSCALE,BZERO");
1359  }
1360  }
1361 
1362  if (scale.bitpix > 0 && !std::numeric_limits<T>::is_integer) {
1363  fits_write_key_lng(fits, "BLANK", scale.blank, "Value for undefined pixels", &status);
1364  fits_write_key_lng(fits, "ZDITHER0", options.scaling.seed, "Dithering seed", &status);
1365  fits_write_key_str(fits, "ZQUANTIZ", "SUBTRACTIVE_DITHER_1", "Dithering algorithm", &status);
1366  if (behavior & AUTO_CHECK) {
1367  LSST_FITS_CHECK_STATUS(*this, "Writing [Z]BLANK");
1368  }
1369  }
1370 
1371  // cfitsio says this is deprecated, but Pan-STARRS found that it was sometimes necessary, writing:
1372  // "This forces a re-scan of the header to ensure everything's kosher.
1373  // Without this, compressed HDUs have been written out with PCOUNT=0 and TFORM1 not correctly set."
1374  fits_set_hdustruc(fits, &status);
1375  if (behavior & AUTO_CHECK) {
1376  LSST_FITS_CHECK_STATUS(*this, "Finalizing header");
1377  }
1378 }
1379 
1380 namespace {
1381 
1383 template <typename T, class Enable = void>
1384 struct NullValue {
1385  static T constexpr value = 0;
1386 };
1387 
1389 template <typename T>
1390 struct NullValue<T, typename std::enable_if<std::numeric_limits<T>::has_quiet_NaN>::type> {
1391  static T constexpr value = std::numeric_limits<T>::quiet_NaN();
1392 };
1393 
1394 } // namespace
1395 
1396 template <typename T>
1397 void Fits::readImageImpl(int nAxis, T *data, long *begin, long *end, long *increment) {
1398  T null = NullValue<T>::value;
1399  int anyNulls = 0;
1400  fits_read_subset(reinterpret_cast<fitsfile *>(fptr), FitsType<T>::CONSTANT, begin, end, increment,
1401  reinterpret_cast<void *>(&null), data, &anyNulls, &status);
1402  if (behavior & AUTO_CHECK) LSST_FITS_CHECK_STATUS(*this, "Reading image");
1403 }
1404 
1406  int nAxis = 0;
1407  fits_get_img_dim(reinterpret_cast<fitsfile *>(fptr), &nAxis, &status);
1408  if (behavior & AUTO_CHECK) LSST_FITS_CHECK_STATUS(*this, "Getting NAXIS");
1409  return nAxis;
1410 }
1411 
1412 void Fits::getImageShapeImpl(int maxDim, long *nAxes) {
1413  fits_get_img_size(reinterpret_cast<fitsfile *>(fptr), maxDim, nAxes, &status);
1414  if (behavior & AUTO_CHECK) LSST_FITS_CHECK_STATUS(*this, "Getting NAXES");
1415 }
1416 
1417 template <typename T>
1419  int imageType = 0;
1420  fits_get_img_equivtype(reinterpret_cast<fitsfile *>(fptr), &imageType, &status);
1421  if (behavior & AUTO_CHECK) LSST_FITS_CHECK_STATUS(*this, "Getting image type");
1423  if (imageType < 0) {
1424  return false; // can't represent floating-point with integer
1425  }
1427  if (isFitsImageTypeSigned(imageType)) {
1428  return FitsBitPix<T>::CONSTANT >= imageType;
1429  } else {
1430  // need extra bits to safely convert unsigned to signed
1431  return FitsBitPix<T>::CONSTANT > imageType;
1432  }
1433  } else {
1434  if (!isFitsImageTypeSigned(imageType)) {
1435  return FitsBitPix<T>::CONSTANT >= imageType;
1436  } else if (imageType == LONGLONG_IMG) {
1437  // workaround for CFITSIO not recognizing uint64 as
1438  // unsigned
1439  return FitsBitPix<T>::CONSTANT >= imageType;
1440  } else {
1441  return false;
1442  }
1443  }
1444  }
1445  // we allow all conversions to float and double, even if they lose precision
1446  return true;
1447 }
1448 
1450  int bitpix = 0;
1451  fits_get_img_equivtype(reinterpret_cast<fitsfile *>(fptr), &bitpix, &status);
1452  if (behavior & AUTO_CHECK) LSST_FITS_CHECK_STATUS(*this, "Getting image type");
1453  // FITS' 'BITPIX' key is the number of bits in a pixel, but negative for
1454  // floats. But the above CFITSIO routine adds support for unsigned
1455  // integers by checking BZERO for an offset as well. So the 'bitpix' value
1456  // we get back from that should be the raw value for signed integers and
1457  // floats, but may be something else (still positive) for unsigned, and
1458  // hence we'll compare to some FITSIO constants to be safe looking at
1459  // integers.
1460  if (bitpix < 0) {
1461  return "float" + std::to_string(-bitpix);
1462  }
1463  switch (bitpix) {
1464  case BYTE_IMG: return "uint8";
1465  case SBYTE_IMG: return "int8";
1466  case SHORT_IMG: return "int16";
1467  case USHORT_IMG: return "uint16";
1468  case LONG_IMG: return "int32";
1469  case ULONG_IMG: return "uint32";
1470  case LONGLONG_IMG: return "int64";
1471  }
1472  throw LSST_EXCEPT(
1473  FitsError,
1474  (boost::format("Unrecognized BITPIX value: %d") % bitpix).str()
1475  );
1476 }
1477 
1479  auto fits = reinterpret_cast<fitsfile *>(fptr);
1480  int compType = 0; // cfitsio compression type
1481  fits_get_compression_type(fits, &compType, &status);
1482  if (behavior & AUTO_CHECK) {
1483  LSST_FITS_CHECK_STATUS(*this, "Getting compression type");
1484  }
1485 
1486  ImageCompressionOptions::Tiles tiles = ndarray::allocate(MAX_COMPRESS_DIM);
1487  fits_get_tile_dim(fits, tiles.getNumElements(), tiles.getData(), &status);
1488  if (behavior & AUTO_CHECK) {
1489  LSST_FITS_CHECK_STATUS(*this, "Getting tile dimensions");
1490  }
1491 
1492  float quantizeLevel;
1493  fits_get_quantize_level(fits, &quantizeLevel, &status);
1494  if (behavior & AUTO_CHECK) {
1495  LSST_FITS_CHECK_STATUS(*this, "Getting quantizeLevel");
1496  }
1497 
1498  return ImageCompressionOptions(compressionAlgorithmFromCfitsio(compType), tiles, quantizeLevel);
1499 }
1500 
1502  auto fits = reinterpret_cast<fitsfile *>(fptr);
1503  fits_unset_compression_request(fits, &status); // wipe the slate and start over
1505  allowImageCompression ? comp.algorithm : ImageCompressionOptions::NONE;
1506  fits_set_compression_type(fits, compressionAlgorithmToCfitsio(algorithm), &status);
1507  if (behavior & AUTO_CHECK) {
1508  LSST_FITS_CHECK_STATUS(*this, "Setting compression type");
1509  }
1510 
1511  if (algorithm == ImageCompressionOptions::NONE) {
1512  // Nothing else worth doing
1513  return;
1514  }
1515 
1516  fits_set_tile_dim(fits, comp.tiles.getNumElements(), comp.tiles.getData(), &status);
1517  if (behavior & AUTO_CHECK) {
1518  LSST_FITS_CHECK_STATUS(*this, "Setting tile dimensions");
1519  }
1520 
1522  fits_set_quantize_level(fits, comp.quantizeLevel, &status);
1523  if (behavior & AUTO_CHECK) {
1524  LSST_FITS_CHECK_STATUS(*this, "Setting quantization level");
1525  }
1526  }
1527 }
1528 
1529 void setAllowImageCompression(bool allow) { allowImageCompression = allow; }
1530 
1531 bool getAllowImageCompression() { return allowImageCompression; }
1532 
1533 // ---- Manipulating files ----------------------------------------------------------------------------------
1534 
1535 Fits::Fits(std::string const &filename, std::string const &mode, int behavior_)
1536  : fptr(0), status(0), behavior(behavior_) {
1537  if (mode == "r" || mode == "rb") {
1538  fits_open_file(reinterpret_cast<fitsfile **>(&fptr), const_cast<char *>(filename.c_str()), READONLY,
1539  &status);
1540  } else if (mode == "w" || mode == "wb") {
1541  boost::filesystem::remove(filename); // cfitsio doesn't like over-writing files
1542  fits_create_file(reinterpret_cast<fitsfile **>(&fptr), const_cast<char *>(filename.c_str()), &status);
1543  } else if (mode == "a" || mode == "ab") {
1544  fits_open_file(reinterpret_cast<fitsfile **>(&fptr), const_cast<char *>(filename.c_str()), READWRITE,
1545  &status);
1546  int nHdu = 0;
1547  fits_get_num_hdus(reinterpret_cast<fitsfile *>(fptr), &nHdu, &status);
1548  fits_movabs_hdu(reinterpret_cast<fitsfile *>(fptr), nHdu, NULL, &status);
1549  if ((behavior & AUTO_CHECK) && (behavior & AUTO_CLOSE) && (status) && (fptr)) {
1550  // We're about to throw an exception, and the destructor won't get called
1551  // because we're in the constructor, so cleanup here first.
1552  int tmpStatus = 0;
1553  fits_close_file(reinterpret_cast<fitsfile *>(fptr), &tmpStatus);
1554  }
1555  } else {
1556  throw LSST_EXCEPT(
1557  FitsError,
1558  (boost::format("Invalid mode '%s' given when opening file '%s'") % mode % filename).str());
1559  }
1560  if (behavior & AUTO_CHECK) {
1561  LSST_FITS_CHECK_STATUS(*this, boost::format("Opening file '%s' with mode '%s'") % filename % mode);
1562  }
1563 }
1564 
1565 Fits::Fits(MemFileManager &manager, std::string const &mode, int behavior_)
1566  : fptr(0), status(0), behavior(behavior_) {
1567  typedef void *(*Reallocator)(void *, std::size_t);
1568  // It's a shame this logic is essentially a duplicate of above, but the innards are different enough
1569  // we can't really reuse it.
1570  if (mode == "r" || mode == "rb") {
1571  fits_open_memfile(reinterpret_cast<fitsfile **>(&fptr), "unused", READONLY, &manager._ptr,
1572  &manager._len, 0, 0, // no reallocator or deltasize necessary for READONLY
1573  &status);
1574  } else if (mode == "w" || mode == "wb") {
1575  Reallocator reallocator = 0;
1576  if (manager._managed) reallocator = &std::realloc;
1577  fits_create_memfile(reinterpret_cast<fitsfile **>(&fptr), &manager._ptr, &manager._len, 0,
1578  reallocator, // use default deltasize
1579  &status);
1580  } else if (mode == "a" || mode == "ab") {
1581  Reallocator reallocator = 0;
1582  if (manager._managed) reallocator = &std::realloc;
1583  fits_open_memfile(reinterpret_cast<fitsfile **>(&fptr), "unused", READWRITE, &manager._ptr,
1584  &manager._len, 0, reallocator, &status);
1585  int nHdu = 0;
1586  fits_get_num_hdus(reinterpret_cast<fitsfile *>(fptr), &nHdu, &status);
1587  fits_movabs_hdu(reinterpret_cast<fitsfile *>(fptr), nHdu, NULL, &status);
1588  if ((behavior & AUTO_CHECK) && (behavior & AUTO_CLOSE) && (status) && (fptr)) {
1589  // We're about to throw an exception, and the destructor won't get called
1590  // because we're in the constructor, so cleanup here first.
1591  int tmpStatus = 0;
1592  fits_close_file(reinterpret_cast<fitsfile *>(fptr), &tmpStatus);
1593  }
1594  } else {
1595  throw LSST_EXCEPT(FitsError,
1596  (boost::format("Invalid mode '%s' given when opening memory file at '%s'") % mode %
1597  manager._ptr)
1598  .str());
1599  }
1600  if (behavior & AUTO_CHECK) {
1602  *this, boost::format("Opening memory file at '%s' with mode '%s'") % manager._ptr % mode);
1603  }
1604 }
1605 
1607  fits_close_file(reinterpret_cast<fitsfile *>(fptr), &status);
1608  fptr = nullptr;
1609 }
1610 
1614  auto combined = std::make_shared<daf::base::PropertyList>();
1615  bool const asScalar = true;
1616  for (auto const &name : first->getOrderedNames()) {
1617  auto const iscv = isCommentIsValid(*first, name);
1618  if (iscv.isComment) {
1619  if (iscv.isValid) {
1620  combined->add<std::string>(name, first->getArray<std::string>(name));
1621  }
1622  } else {
1623  combined->copy(name, first, name, asScalar);
1624  }
1625  }
1626  for (auto const &name : second->getOrderedNames()) {
1627  auto const iscv = isCommentIsValid(*second, name);
1628  if (iscv.isComment) {
1629  if (iscv.isValid) {
1630  combined->add<std::string>(name, second->getArray<std::string>(name));
1631  }
1632  } else {
1633  // `copy` will replace an item, even if has a different type, so no need to call `remove`
1634  combined->copy(name, second, name, asScalar);
1635  }
1636  }
1637  return combined;
1638 }
1639 
1642  fp.setHdu(hdu);
1643  return readMetadata(fp, strip);
1644 }
1645 
1648  fp.setHdu(hdu);
1649  return readMetadata(fp, strip);
1650 }
1651 
1653  auto metadata = std::make_shared<lsst::daf::base::PropertyList>();
1654  fitsfile.readMetadata(*metadata, strip);
1655  // if INHERIT=T, we want to also include header entries from the primary HDU
1656  int oldHdu = fitsfile.getHdu();
1657  if (oldHdu != 0 && metadata->exists("INHERIT")) {
1658  bool inherit = false;
1659  if (metadata->typeOf("INHERIT") == typeid(std::string)) {
1660  inherit = (metadata->get<std::string>("INHERIT") == "T");
1661  } else {
1662  inherit = metadata->get<bool>("INHERIT");
1663  }
1664  if (strip) metadata->remove("INHERIT");
1665  if (inherit) {
1666  HduMoveGuard guard(fitsfile, 0);
1667  // Combine the metadata from the primary HDU with the metadata from the specified HDU,
1668  // with non-comment values from the specified HDU superseding those in the primary HDU
1669  // and comments from the specified HDU appended to comments from the primary HDU
1670  auto primaryHduMetadata = std::make_shared<daf::base::PropertyList>();
1671  fitsfile.readMetadata(*primaryHduMetadata, strip);
1672  metadata = combineMetadata(primaryHduMetadata, metadata);
1673  } else {
1674  // Purge invalid values
1675  auto const emptyMetadata = std::make_shared<lsst::daf::base::PropertyList>();
1676  metadata = combineMetadata(metadata, emptyMetadata);
1677  }
1678  }
1679  return metadata;
1680 }
1681 
1682 
1683 HduMoveGuard::HduMoveGuard(Fits & fits, int hdu, bool relative) :
1684  _fits(fits),
1685  _oldHdu(_fits.getHdu()),
1686  _enabled(true)
1687 {
1688  _fits.setHdu(hdu, relative);
1689 }
1690 
1692  if (!_enabled) {
1693  return;
1694  }
1695  int status = 0;
1696  std::swap(status, _fits.status); // unset error indicator, but remember the old status
1697  try {
1698  _fits.setHdu(_oldHdu);
1699  } catch (...) {
1700  LOGL_WARN(
1701  "afw.fits",
1702  makeErrorMessage(_fits.fptr, _fits.status, "Failed to move back to HDU %d").c_str(),
1703  _oldHdu
1704  );
1705  }
1706  std::swap(status, _fits.status); // reset the old status
1707 }
1708 
1710  auto fits = reinterpret_cast<fitsfile *>(fptr);
1711  if (getHdu() != 0 || countHdus() == 1) {
1712  return false; // Can't possibly be the PHU leading a compressed image
1713  }
1714  // Check NAXIS = 0
1715  int naxis;
1716  fits_get_img_dim(fits, &naxis, &status);
1717  if (behavior & AUTO_CHECK) {
1718  LSST_FITS_CHECK_STATUS(*this, "Checking NAXIS of PHU");
1719  }
1720  if (naxis != 0) {
1721  return false;
1722  }
1723  // Check first extension (and move back there when we're done if we're not compressed)
1724  HduMoveGuard move(*this, 1);
1725  bool isCompressed = fits_is_compressed_image(fits, &status);
1726  if (behavior & AUTO_CHECK) {
1727  LSST_FITS_CHECK_STATUS(*this, "Checking compression");
1728  }
1729  if (isCompressed) {
1730  move.disable();
1731  }
1732  return isCompressed;
1733 }
1734 
1736  : compression(fits::compressionAlgorithmFromString(config.get<std::string>("compression.algorithm")),
1737  std::vector<long>{config.getAsInt64("compression.columns"),
1738  config.getAsInt64("compression.rows")},
1739  config.getAsDouble("compression.quantizeLevel")),
1740  scaling(fits::scalingAlgorithmFromString(config.get<std::string>("scaling.algorithm")),
1741  config.getAsInt("scaling.bitpix"),
1742  config.exists("scaling.maskPlanes") ? config.getArray<std::string>("scaling.maskPlanes")
1744  config.getAsInt("scaling.seed"), config.getAsDouble("scaling.quantizeLevel"),
1745  config.getAsDouble("scaling.quantizePad"), config.get<bool>("scaling.fuzz"),
1746  config.getAsDouble("scaling.bscale"), config.getAsDouble("scaling.bzero")) {}
1747 
1748 namespace {
1749 
1750 template <typename T>
1751 void validateEntry(daf::base::PropertySet &output, daf::base::PropertySet const &input,
1752  std::string const &name, T defaultValue) {
1753  output.add(name, input.get<T>(name, defaultValue));
1754 }
1755 
1756 template <typename T>
1757 void validateEntry(daf::base::PropertySet &output, daf::base::PropertySet const &input,
1758  std::string const &name, std::vector<T> defaultValue) {
1759  output.add(name, input.exists(name) ? input.getArray<T>(name) : defaultValue);
1760 }
1761 
1762 } // anonymous namespace
1763 
1765  auto validated = std::make_shared<daf::base::PropertySet>();
1766 
1767  validateEntry(*validated, config, "compression.algorithm", std::string("NONE"));
1768  validateEntry(*validated, config, "compression.columns", 0);
1769  validateEntry(*validated, config, "compression.rows", 1);
1770  validateEntry(*validated, config, "compression.quantizeLevel", 0.0);
1771 
1772  validateEntry(*validated, config, "scaling.algorithm", std::string("NONE"));
1773  validateEntry(*validated, config, "scaling.bitpix", 0);
1774  validateEntry(*validated, config, "scaling.maskPlanes", std::vector<std::string>{"NO_DATA"});
1775  validateEntry(*validated, config, "scaling.seed", 1);
1776  validateEntry(*validated, config, "scaling.quantizeLevel", 5.0);
1777  validateEntry(*validated, config, "scaling.quantizePad", 10.0);
1778  validateEntry(*validated, config, "scaling.fuzz", true);
1779  validateEntry(*validated, config, "scaling.bscale", 1.0);
1780  validateEntry(*validated, config, "scaling.bzero", 0.0);
1781 
1782  // Check for additional entries that we don't support (e.g., from typos)
1783  for (auto const &name : config.names(false)) {
1784  if (!validated->exists(name)) {
1786  os << "Invalid image write option: " << name;
1788  }
1789  }
1790 
1791  return validated;
1792 }
1793 
1794 #define INSTANTIATE_KEY_OPS(r, data, T) \
1795  template void Fits::updateKey(std::string const &, T const &, std::string const &); \
1796  template void Fits::writeKey(std::string const &, T const &, std::string const &); \
1797  template void Fits::updateKey(std::string const &, T const &); \
1798  template void Fits::writeKey(std::string const &, T const &); \
1799  template void Fits::updateColumnKey(std::string const &, int, T const &, std::string const &); \
1800  template void Fits::writeColumnKey(std::string const &, int, T const &, std::string const &); \
1801  template void Fits::updateColumnKey(std::string const &, int, T const &); \
1802  template void Fits::writeColumnKey(std::string const &, int, T const &); \
1803  template void Fits::readKey(std::string const &, T &);
1804 
1805 #define INSTANTIATE_IMAGE_OPS(r, data, T) \
1806  template void Fits::writeImageImpl(T const *, int); \
1807  template void Fits::writeImage(image::ImageBase<T> const &, ImageWriteOptions const &, \
1808  std::shared_ptr<daf::base::PropertySet const>, \
1809  std::shared_ptr<image::Mask<image::MaskPixel> const>); \
1810  template void Fits::readImageImpl(int, T *, long *, long *, long *); \
1811  template bool Fits::checkImageType<T>(); \
1812  template int getBitPix<T>();
1813 
1814 #define INSTANTIATE_TABLE_OPS(r, data, T) \
1815  template int Fits::addColumn<T>(std::string const &ttype, int size); \
1816  template int Fits::addColumn<T>(std::string const &ttype, int size, std::string const &comment);
1817 #define INSTANTIATE_TABLE_ARRAY_OPS(r, data, T) \
1818  template void Fits::writeTableArray(std::size_t row, int col, int nElements, T const *value); \
1819  template void Fits::readTableArray(std::size_t row, int col, int nElements, T *value);
1820 
1821 // ----------------------------------------------------------------------------------------------------------
1822 // ---- Explicit instantiation ------------------------------------------------------------------------------
1823 // ----------------------------------------------------------------------------------------------------------
1824 
1825 #define KEY_TYPES \
1826  (bool)(unsigned char)(short)(unsigned short)(int)(unsigned int)(long)(unsigned long)(LONGLONG)( \
1827  float)(double)(std::complex<float>)(std::complex<double>)(std::string)
1828 
1829 #define COLUMN_TYPES \
1830  (bool)(std::string)(std::int8_t)(std::uint8_t)(std::int16_t)(std::uint16_t)(std::int32_t)(std::uint32_t) \
1831  (std::int64_t)(float)(double)(lsst::geom::Angle)(std::complex<float>)(std::complex<double>)
1832 
1833 #define COLUMN_ARRAY_TYPES \
1834  (bool)(char)(std::uint8_t)(std::int16_t)(std::uint16_t)(std::int32_t)(std::uint32_t)(std::int64_t)( \
1835  float)(double)(lsst::geom::Angle)(std::complex<float>)(std::complex<double>)
1836 
1837 #define IMAGE_TYPES \
1838  (unsigned char)(short)(unsigned short)(int)(unsigned int)(std::int64_t)(std::uint64_t)(float)(double)
1839 
1844 } // namespace fits
1845 } // namespace afw
1846 } // namespace lsst
#define LOGLS_WARN(logger, message)
Log a warn-level message using an iostream-based interface.
Definition: Log.h:648
int col
Definition: CR.cc:144
def write(self, patchRef, catalog)
Write the output.
int getImageDim()
Return the number of dimensions in the current HDU.
Definition: fits.cc:1405
std::size_t countRows()
Return the number of row in a table.
Definition: fits.cc:1142
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
#define INSTANTIATE_TABLE_OPS(r, data, T)
Definition: fits.cc:1814
std::shared_ptr< daf::base::PropertyList > combineMetadata(std::shared_ptr< const daf::base::PropertyList > first, std::shared_ptr< const daf::base::PropertyList > second)
Combine two sets of metadata in a FITS-appropriate fashion.
Definition: fits.cc:1611
std::string getFileName() const
Return the file name associated with the FITS object or "<unknown>" if there is none.
Definition: fits.cc:497
void forEachKey(HeaderIterationFunctor &functor)
Call a polymorphic functor for every key in the header.
Definition: fits.cc:780
T empty(T... args)
ImageCompressionOptions getImageCompression()
Return the current image compression settings.
Definition: fits.cc:1478
ImageCompressionOptions compression
Options controlling compression.
Definition: fits.h:220
T find_first_not_of(T... args)
afw::table::Key< afw::table::Array< MaskPixelT > > mask
void writeMetadata(daf::base::PropertySet const &metadata)
Read a FITS header into a PropertySet or PropertyList.
Definition: fits.cc:1078
#define LSST_FITS_CHECK_STATUS(fitsObj,...)
Throw a FitsError exception if the status of the given Fits object is nonzero.
Definition: fits.h:111
std::vector< T > getArray(std::string const &name) const
Get the vector of values for a property name (possibly hierarchical).
int compressionAlgorithmToCfitsio(ImageCompressionOptions::CompressionAlgorithm algorithm)
Convert ImageCompressionOptions::CompressionAlgorithm to cfitsio.
T swap(T... args)
T front(T... args)
Class for storing ordered metadata with comments.
Definition: PropertyList.h:68
T get(std::string const &name) const
Get the last value for a property name (possibly hierarchical).
long blank
Value for integer images indicating non-finite values.
void readKey(std::string const &key, T &value)
Read a FITS header key into the given reference.
Definition: fits.cc:773
int bitpix
Bits per pixel; negative means floating-point: 8,16,32,64,-32,-64.
void updateColumnKey(std::string const &prefix, int n, T const &value, std::string const &comment)
Update a key of the form XXXXXnnn, where XXXXX is the prefix and nnn is a column number.
Definition: fits.cc:690
Options for writing an image to FITS.
Definition: fits.h:219
std::string const & getComment(std::string const &name) const
Get the comment for a string property name (possibly hierarchical).
Definition: PropertyList.cc:77
def scale(algorithm, min, max=None, frame=None)
Definition: ds9.py:109
T to_string(T... args)
int getBitPix()
Return the cfitsio integer BITPIX code for the given data type.
Definition: fits.cc:489
void readTableScalar(std::size_t row, int col, T &value)
Read an array scalar from a binary table.
Definition: fits.h:595
bool checkImageType()
Return true if the current HDU is compatible with the given pixel type.
Definition: fits.cc:1418
STL namespace.
ImageT val
Definition: CR.cc:146
STL class.
void writeColumnKey(std::string const &prefix, int n, T const &value, std::string const &comment)
Write a key of the form XXXXXnnn, where XXXXX is the prefix and nnn is a column number.
Definition: fits.cc:698
T end(T... args)
#define INSTANTIATE_TABLE_ARRAY_OPS(r, data, T)
Definition: fits.cc:1817
std::vector< std::string > paramNames(bool topLevelOnly=true) const
A variant of names that excludes the names of subproperties.
CompressionAlgorithm
Compression algorithms.
T free(T... args)
daf::base::PropertySet * set
Definition: fits.cc:902
A simple struct that combines the two arguments that must be passed to most cfitsio routines and cont...
Definition: fits.h:297
The base class for all image classed (Image, Mask, MaskedImage, ...)
Definition: ImageBase.h:102
long getTableArraySize(int col)
Return the size of an array column.
Definition: fits.cc:1200
FITS BITPIX header value by C++ type.
void closeFile()
Close a FITS file.
Definition: fits.cc:1606
A class representing an angle.
Definition: Angle.h:127
int seed
Seed for random number generator when fuzzing.
Fits * fits
Definition: FitsWriter.cc:90
void createEmpty()
Create an empty image HDU with NAXIS=0 at the end of the file.
Definition: fits.cc:1223
std::string getImageDType()
Return the numpy dtype equivalent of the image pixel type (e.g.
Definition: fits.cc:1449
void updateKey(std::string const &key, T const &value, std::string const &comment)
Set a FITS header key, editing if it already exists and appending it if not.
Definition: fits.cc:658
STL class.
LSST DM logging module built on log4cxx.
#define KEY_TYPES
Definition: fits.cc:1825
ndarray::Array< T const, N, N > const makeContiguousArray(ndarray::Array< T, N, C > const &array)
Construct a contiguous ndarray.
Definition: fits.h:208
An exception thrown when problems are found when reading or writing FITS files.
Definition: fits.h:36
Tiles tiles
Tile size; a dimension with 0 means infinite (e.g., to specify one row: 0,1)
T push_back(T... args)
Options for tile compression of image pixels.
A base class for image defects.
char * data
Definition: BaseRecord.cc:62
table::Key< int > type
Definition: Detector.cc:163
lsst::geom::Box2I getBBox(ImageOrigin origin=PARENT) const
Definition: ImageBase.h:462
Represent a 2-dimensional array of bitmask pixels.
Definition: Mask.h:77
int64_t getAsInt64(std::string const &name) const
Get the last value for a bool/char/short/int/int64_t property name (possibly hierarchical).
Lifetime-management for memory that goes into FITS memory files.
Definition: fits.h:121
double bzero
Zero-point to apply when reading from FITS.
std::size_t addRows(std::size_t nRows)
Append rows to a table, and return the index of the first new row.
Definition: fits.cc:1132
bool isValid
Definition: fits.cc:398
std::shared_ptr< daf::base::PropertyList > createTrivialWcsMetadata(std::string const &wcsName, lsst::geom::Point2I const &xy0)
Definition: wcsUtils.cc:51
ImageCompressionOptions::CompressionAlgorithm compressionAlgorithmFromString(std::string const &name)
Interpret compression algorithm expressed in string.
int getArea() const
Definition: Box.h:189
Fits()
Default constructor; set all data members to 0.
Definition: fits.h:609
T erase(T... args)
Key< U > key
Definition: Schema.cc:281
T str(T... args)
T isfinite(T... args)
bool exists(std::string const &name) const
Determine if a name (possibly hierarchical) exists.
#define IMAGE_TYPES
Definition: fits.cc:1837
ndarray::Array< long, 1, 1 > Tiles
void createTable()
Create a new binary table extension.
Definition: fits.cc:1100
#define LOGL_WARN(logger, message...)
Log a warn-level message using a varargs/printf style interface.
Definition: Log.h:536
double bscale
Scale to apply when reading from FITS.
#define COLUMN_TYPES
Definition: fits.cc:1829
T infinity(T... args)
void setHdu(int hdu, bool relative=false)
Set the current HDU.
Definition: fits.cc:512
T count(T... args)
T find_last_not_of(T... args)
#define INSTANTIATE_IMAGE_OPS(r, data, T)
Definition: fits.cc:1805
BOOST_PP_SEQ_FOR_EACH(INSTANTIATE_COLUMNVIEW_SCALAR, _, BOOST_PP_TUPLE_TO_SEQ(AFW_TABLE_SCALAR_FIELD_TYPE_N, AFW_TABLE_SCALAR_FIELD_TYPE_TUPLE)) BOOST_PP_SEQ_FOR_EACH(INSTANTIATE_COLUMNVIEW_ARRAY
lsst::geom::Point2I getXY0() const
Return the image&#39;s origin.
Definition: ImageBase.h:340
T find(T... args)
T size(T... args)
STL class.
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
T realloc(T... args)
T name(T... args)
#define INSTANTIATE_KEY_OPS(r, data, T)
Definition: fits.cc:1794
T malloc(T... args)
ImageScalingOptions scaling
Options controlling scaling.
Definition: fits.h:221
static std::shared_ptr< daf::base::PropertySet > validate(daf::base::PropertySet const &config)
Validate a PropertySet.
Definition: fits.cc:1764
void disable()
Disable the guard, leaving the HDU at its current state at destruction.
Definition: fits.h:751
bool checkCompressedImagePhu()
Go to the first image header in the FITS file.
Definition: fits.cc:1709
void setAllowImageCompression(bool allow)
Definition: fits.cc:1529
Class for storing generic metadata.
Definition: PropertySet.h:67
T c_str(T... args)
float quantizeLevel
quantization level: 0.0 = none requires use of GZIP or GZIP_SHUFFLE
void writeTableArray(std::size_t row, int col, int nElements, T const *value)
Write an array value to a binary table.
Definition: fits.cc:1152
void reset()
Return the manager to the same state it would be if default-constructed.
Definition: fits.cc:474
CompressionAlgorithm algorithm
Compresion algorithm to use.
void readTableArray(std::size_t row, int col, int nElements, T *value)
Read an array value from a binary table.
Definition: fits.cc:1175
std::string prefix
Definition: SchemaMapper.cc:79
T isnan(T... args)
T substr(T... args)
ImageCompressionOptions::CompressionAlgorithm compressionAlgorithmFromCfitsio(int cfitsio)
Convert compression algorithm from cfitsio to ImageCompressionOptions::CompressionAlgorithm.
ImageScale determine(image::ImageBase< T > const &image, std::shared_ptr< image::Mask< image::MaskPixel > const > mask=nullptr) const
Determine the scaling for a particular image.
void setImageCompression(ImageCompressionOptions const &options)
Set compression options for writing FITS images.
Definition: fits.cc:1501
#define COLUMN_ARRAY_TYPES
Definition: fits.cc:1833
T quiet_NaN(T... args)
void writeKey(std::string const &key, T const &value, std::string const &comment)
Add a FITS header key to the bottom of the header.
Definition: fits.cc:666
bool fuzz
Fuzz the values when quantising floating-point values?
std::vector< std::string > names(bool topLevelOnly=true) const
Get the names in the PropertySet, optionally including those in subproperties.
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects...
bool getAllowImageCompression()
Definition: fits.cc:1531
int getHdu()
Return the current HDU (0-indexed; 0 is the Primary HDU).
Definition: fits.cc:506
void add(std::string const &name, T const &value)
Append a single value to the vector of values for a property name (possibly hierarchical).
void writeImage(ndarray::Array< T const, N, C > const &array)
Write an ndarray::Array to a FITS image HDU.
Definition: fits.h:477
std::ostream * os
Definition: Schema.cc:746
def writeMetadata(self, dataRefList)
No metadata to write, and not sure how to write it for a list of dataRefs.
std::vector< std::string > getOrderedNames() const
Get the list of property names, in the order they were added.
Definition: PropertyList.cc:81
int addColumn(std::string const &ttype, int size, std::string const &comment)
Add a column to a table.
Definition: fits.cc:1123
daf::base::PropertyList * list
Definition: fits.cc:903
std::string makeErrorMessage(std::string const &fileName="", int status=0, std::string const &msg="")
Return an error message reflecting FITS I/O errors.
Definition: fits.cc:425
py::object result
Definition: _schema.cc:429
RAII scoped guard for moving the HDU in a Fits object.
Definition: fits.h:724
ImageScalingOptions::ScalingAlgorithm scalingAlgorithmFromString(std::string const &name)
Interpret scaling algorithm expressed in string.
std::string const wcsNameForXY0
Definition: ImageBase.h:70
Base class for polymorphic functors used to iterator over FITS key headers.
Definition: fits.h:49
int end
bool strip
Definition: fits.cc:901
void writeTableScalar(std::size_t row, int col, T value)
Write a scalar value to a binary table.
Definition: fits.h:583
void readMetadata(daf::base::PropertySet &metadata, bool strip=false)
Read a FITS header into a PropertySet or PropertyList.
Definition: fits.cc:1070
std::string makeLimitedFitsHeader(lsst::daf::base::PropertySet const &metadata, std::set< std::string > const &excludeNames={})
Format a PropertySet into an FITS header string in a simplistic fashion.
Definition: fits.cc:456
T compare(T... args)
bool isComment
Definition: fits.cc:397
T reserve(T... args)
Reports errors that are due to events beyond the control of the program.
Definition: Runtime.h:104
const int DEFAULT_HDU
Specify that the default HDU should be read.
Definition: fitsDefaults.h:18
std::shared_ptr< detail::PixelArrayBase > toFits(ndarray::Array< T const, 2, 2 > const &image, bool forceNonfiniteRemoval, bool fuzz=true, ndarray::Array< long, 1 > const &tiles=ndarray::Array< long, 1, 1 >(), int seed=1) const
Convert to an array of pixel values to write to FITS.
int row
Definition: CR.cc:145
ImageWriteOptions(image::Image< T > const &image)
Construct with default options for images.
Definition: fits.h:225
int countHdus()
Return the number of HDUs in the file.
Definition: fits.cc:533
Scale to apply to image.