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