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