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