11 #include "boost/regex.hpp"
12 #include "boost/multi_index_container.hpp"
13 #include "boost/multi_index/sequenced_index.hpp"
14 #include "boost/multi_index/ordered_index.hpp"
15 #include "boost/multi_index/hashed_index.hpp"
16 #include "boost/multi_index/member.hpp"
35 template <std::
string FitsSchemaItem::*Member>
36 struct SetFitsSchemaString {
37 void operator()(FitsSchemaItem &item) { item.*Member = _v; }
38 explicit SetFitsSchemaString(
std::string const &v) : _v(v) {}
52 typedef boost::multi_index_container<
54 boost::multi_index::indexed_by<
55 boost::multi_index::ordered_non_unique<
56 boost::multi_index::member<FitsSchemaItem, int, &FitsSchemaItem::column>>,
57 boost::multi_index::ordered_non_unique<
58 boost::multi_index::member<FitsSchemaItem, int, &FitsSchemaItem::bit>>,
59 boost::multi_index::hashed_unique<
60 boost::multi_index::member<FitsSchemaItem, std::string, &FitsSchemaItem::ttype>>,
61 boost::multi_index::sequenced<>>>
65 typedef SetFitsSchemaString<&FitsSchemaItem::ttype>
SetTTYPE;
66 typedef SetFitsSchemaString<&FitsSchemaItem::tform>
SetTFORM;
67 typedef SetFitsSchemaString<&FitsSchemaItem::tccls>
SetTCCLS;
68 typedef SetFitsSchemaString<&FitsSchemaItem::tunit>
SetTUNIT;
69 typedef SetFitsSchemaString<&FitsSchemaItem::doc>
SetDoc;
101 : _impl(
std::make_shared<
Impl>()) {
104 if (!metadata.
exists(
"AFW_TYPE")) {
107 _impl->version = metadata.
get(
"AFW_TABLE_VERSION", _impl->version);
108 _impl->type = metadata.
get(
"AFW_TYPE", _impl->type);
110 metadata.
remove(
"AFW_TABLE_VERSION");
113 metadata.
remove(
"AFW_TYPE");
117 _impl->archiveHdu = metadata.
get(
"AR_HDU", -1);
118 if (_impl->archiveHdu > 0) {
121 metadata.
remove(
"AR_HDU");
130 if (pos == std::string::npos) {
132 (
boost::format(
"Malformed alias definition: '%s'") % (*i)).str());
134 _impl->schema.getAliasMap()->set(i->substr(0, pos), i->substr(pos + 1, std::string::npos));
143 if (_impl->version == 0) {
155 _impl->schema.getAliasMap()->set(oldSlotKeys[i].
second,
target);
158 metadata.
remove(oldSlotKeys[i].
first +
"_ERR_SLOT");
159 metadata.
remove(oldSlotKeys[i].
first +
"_FLAG_SLOT");
168 if (
key->compare(0, 5,
"TTYPE") == 0) {
170 auto iter = _impl->byColumn().lower_bound(column);
171 if (
iter == _impl->byColumn().end() ||
iter->column != column) {
176 if (
iter->doc.empty()) {
182 }
else if (
key->compare(0, 5,
"TFLAG") == 0) {
184 auto iter = _impl->byBit().lower_bound(bit);
185 if (
iter == _impl->byBit().end() ||
iter->bit != bit) {
190 if (
iter->doc.empty()) {
196 }
else if (
key->compare(0, 4,
"TDOC") == 0) {
198 auto iter = _impl->byColumn().lower_bound(column);
199 if (
iter == _impl->byColumn().end() ||
iter->column != column) {
206 }
else if (
key->compare(0, 5,
"TFDOC") == 0) {
208 auto iter = _impl->byBit().lower_bound(bit);
209 if (
iter == _impl->byBit().end() ||
iter->bit != bit) {
216 }
else if (
key->compare(0, 5,
"TUNIT") == 0) {
218 auto iter = _impl->byColumn().lower_bound(column);
219 if (
iter == _impl->byColumn().end() ||
iter->column != column) {
226 }
else if (
key->compare(0, 5,
"TCCLS") == 0) {
228 auto iter = _impl->byColumn().lower_bound(column);
229 if (
iter == _impl->byColumn().end() ||
iter->column != column) {
236 }
else if (
key->compare(0, 5,
"TFORM") == 0) {
238 auto iter = _impl->byColumn().lower_bound(column);
239 if (
iter == _impl->byColumn().end() ||
iter->column != column) {
246 }
else if (
key->compare(0, 5,
"TZERO") == 0) {
250 }
else if (
key->compare(0, 5,
"TSCAL") == 0) {
254 }
else if (
key->compare(0, 5,
"TNULL") == 0) {
258 }
else if (
key->compare(0, 5,
"TDISP") == 0) {
266 _impl->flagColumn = metadata.
get(
"FLAGCOL", 0);
267 if (_impl->flagColumn > 0) {
269 metadata.
remove(
"FLAGCOL");
272 auto iter = _impl->byColumn().find(_impl->flagColumn);
273 if (
iter == _impl->byColumn().end()) {
276 (
boost::format(
"Column for flag data not found; FLAGCOL=%d") % _impl->flagColumn).str());
281 static boost::regex
const regex(
"(\\d+)?X\\(?(\\d)*\\)?", boost::regex::perl);
283 if (!boost::regex_match(
iter->tform,
m, regex)) {
286 (
boost::format(
"Invalid TFORM key for flags column: '%s'") %
iter->tform).str());
292 _impl->flagKeys.resize(
nFlags);
293 _impl->flagWorkspace.reset(
new bool[
nFlags]);
296 _impl->byColumn().erase(
iter);
309 int oldHdu =
fits.getHdu();
310 if (_impl->archiveHdu < 0) _impl->archiveHdu = oldHdu + 1;
312 fits.setHdu(_impl->archiveHdu);
319 _impl->archiveHdu = -1;
327 auto iter = _impl->byName().find(ttype);
328 if (
iter == _impl->byName().end()) {
335 auto iter = _impl->byColumn().lower_bound(column);
336 if (
iter == _impl->byColumn().end() ||
iter->column != column) {
343 auto iter = _impl->byColumn().lower_bound(item->
column);
344 assert(
iter != _impl->byColumn().end() &&
iter->column == item->
column);
345 _impl->byColumn().erase(
iter);
349 auto iter = _impl->byName().find(ttype);
350 if (
iter != _impl->byName().end() &&
iter->ttype == ttype) {
351 _impl->byName().erase(
iter);
356 auto iter = _impl->byColumn().lower_bound(column);
357 if (
iter != _impl->byColumn().end() &&
iter->column == column) {
358 _impl->byColumn().erase(
iter);
365 _impl->readers.push_back(
std::move(reader));
370 template <
typename T>
379 : _column(item.column), _key(
schema.addField<T>(item.ttype, item.doc, item.tunit, base)),
380 _cache(), _cacheFirstRow(0)
387 if (_key.getElementCount() == 1u) {
388 std::size_t nElements = nRows*_key.getElementCount();
389 _cache.resize(nElements);
390 _cacheFirstRow = firstRow;
391 fits.readTableArray(firstRow, _column, nElements, &_cache.front());
397 if (_cache.empty()) {
398 fits.readTableArray(
row, _column, _key.getElementCount(), record.getElement(_key));
400 assert(
row >= _cacheFirstRow);
402 assert(offset < _cache.size());
403 std::copy_n(_cache.begin() + offset, _key.getElementCount(), record.getElement(_key));
415 class AngleReader :
public FitsColumnReader {
418 Schema &
schema, FitsSchemaItem
const &item,
419 FieldBase<lsst::geom::Angle>
const &base = FieldBase<lsst::geom::Angle>()) {
423 AngleReader(Schema &
schema, FitsSchemaItem
const &item, FieldBase<lsst::geom::Angle>
const &base)
424 : _column(item.column), _key(
schema.addField<
lsst::
geom::
Angle>(item.ttype, item.doc,
"", base)) {
429 if (!item.tunit.empty() && item.tunit !=
"rad") {
431 "Angle fields must be persisted in radians (TUNIT='rad').");
436 assert(_key.getElementCount() == 1u);
437 _cache.resize(nRows);
438 _cacheFirstRow = firstRow;
439 fits.readTableArray(firstRow, _column, nRows, &_cache.front());
444 if (_cache.empty()) {
446 fits.readTableScalar(
row, _column, tmp);
449 assert(
row >= _cacheFirstRow);
451 assert(offset < _cache.size());
458 Key<lsst::geom::Angle> _key;
463 class StringReader :
public FitsColumnReader {
469 StringReader(Schema &
schema, FitsSchemaItem
const &item,
int size)
470 : _column(item.column),
471 _key(
schema.addField<
std::string>(item.ttype, item.doc, item.tunit, size)),
472 _isVariableLength(size == 0) {}
477 fits.readTableScalar(
row, _column, s, _isVariableLength);
483 Key<std::string> _key;
484 bool _isVariableLength;
487 template <
typename T>
488 class VariableLengthArrayReader :
public FitsColumnReader {
494 VariableLengthArrayReader(Schema &
schema, FitsSchemaItem
const &item)
495 : _column(item.column), _key(
schema.addField<Array<T>>(item.ttype, item.doc, item.tunit, 0)) {}
499 int size =
fits.getTableArraySize(
row, _column);
500 ndarray::Array<T, 1, 1> array = ndarray::allocate(size);
501 fits.readTableArray(
row, _column, size, array.getData());
502 record.set(_key, array);
512 template <
typename T>
513 class PointConversionReader :
public FitsColumnReader {
519 PointConversionReader(Schema &
schema, FitsSchemaItem
const &item)
520 : _column(item.column), _key(PointKey<T>::addFields(
schema, item.ttype, item.doc, item.tunit)) {}
525 fits.readTableArray(
row, _column, 2, buffer.
data());
536 class CoordConversionReader :
public FitsColumnReader {
542 CoordConversionReader(Schema &
schema, FitsSchemaItem
const &item)
543 : _column(item.column), _key(CoordKey::addFields(
schema, item.ttype, item.doc)) {}
548 fits.readTableArray(
row, _column, 2, buffer.
data());
559 class MomentsConversionReader :
public FitsColumnReader {
565 MomentsConversionReader(Schema &
schema, FitsSchemaItem
const &item)
566 : _column(item.column),
572 fits.readTableArray(
row, _column, 3, buffer.
data());
573 record.set(_key, geom::ellipses::Quadrupole(buffer[0], buffer[1], buffer[2],
false));
584 template <
typename T,
int N>
585 class CovarianceConversionReader :
public FitsColumnReader {
588 static boost::regex
const regex(
"(.*)(\\^(\\d+))?", boost::regex::perl);
590 if (!boost::regex_match(oldUnits,
m, regex)) {
603 CovarianceConversionReader(Schema &
schema, FitsSchemaItem
const &item,
605 : _column(item.column),
607 _key(CovarianceMatrixKey<T, N>::addFields(
schema, item.ttype, names, guessUnits(item.tunit))),
613 for (
int i = 0; i <
_size; ++i) {
614 for (
int j = i; j <
_size; ++j) {
623 CovarianceMatrixKey<T, N> _key;
631 static boost::regex
const regex(
"(\\d+)?([PQ])?(\\u)\\(?(\\d)*\\)?", boost::regex::perl);
634 if (!boost::regex_match(item.tform,
m, regex)) {
641 char code =
m[3].str()[0];
651 if (item.tccls ==
"Array") {
652 return StandardReader<Array<std::uint8_t>>::make(
schema, item, size);
654 return StandardReader<std::uint8_t>::make(
schema, item);
657 return VariableLengthArrayReader<std::uint8_t>::make(
schema, item);
659 return StandardReader<Array<std::uint8_t>>::make(
schema, item, size);
664 if (item.tccls ==
"Array") {
665 return StandardReader<Array<std::uint16_t>>::make(
schema, item, size);
667 return StandardReader<std::uint16_t>::make(
schema, item);
670 return VariableLengthArrayReader<std::uint16_t>::make(
schema, item);
672 return StandardReader<Array<std::uint16_t>>::make(
schema, item, size);
675 return VariableLengthArrayReader<std::int32_t>::make(
schema, item);
677 if (item.tccls ==
"Point") {
678 return PointConversionReader<std::int32_t>::make(
schema, item);
680 if (size > 1 || item.tccls ==
"Array") {
681 return StandardReader<Array<std::int32_t>>::make(
schema, item, size);
683 return StandardReader<std::int32_t>::make(
schema, item);
686 return StandardReader<std::int64_t>::make(
schema, item);
690 return VariableLengthArrayReader<float>::make(
schema, item);
693 if (item.tccls ==
"Array") {
694 return StandardReader<Array<float>>::make(
schema, item, 1);
698 return StandardReader<float>::make(
schema, item);
700 if (size == 3 && item.tccls ==
"Covariance(Point)") {
702 return CovarianceConversionReader<float, 2>::make(
schema, item, names);
704 if (size == 6 && item.tccls ==
"Covariance(Moments)") {
706 return CovarianceConversionReader<float, 3>::make(
schema, item, names);
708 if (item.tccls ==
"Covariance") {
709 double v = 0.5 * (
std::sqrt(1 + 8 * size) - 1);
711 if (n * (n + 1) != size * 2) {
712 throw LSST_EXCEPT(afw::fits::FitsError,
"Covariance field has invalid size.");
715 for (
int i = 0; i < n; ++i) {
718 return CovarianceConversionReader<float, Eigen::Dynamic>::make(
schema, item, names);
720 return StandardReader<Array<float>>::make(
schema, item, size);
723 return VariableLengthArrayReader<double>::make(
schema, item);
726 if (item.tccls ==
"Angle") {
727 return AngleReader::make(
schema, item);
729 if (item.tccls ==
"Array") {
730 return StandardReader<Array<double>>::make(
schema, item, 1);
732 return StandardReader<double>::make(
schema, item);
735 if (item.tccls ==
"Point") {
736 return PointConversionReader<double>::make(
schema, item);
738 if (item.tccls ==
"Coord") {
739 return CoordConversionReader::make(
schema, item);
742 if (size == 3 && item.tccls ==
"Moments") {
743 return MomentsConversionReader::make(
schema, item);
745 return StandardReader<Array<double>>::make(
schema, item, size);
748 return StringReader::make(
schema, item, size);
758 bool isInstFlux(FitsSchemaItem
const & item) {
763 if (!
includes(item.ttype,
"flux"))
return false;
764 if (
includes(item.ttype,
"modelfit_CModel") && item.tunit.empty()) {
771 std::transform(units.begin(), units.end(), units.begin(), [](
char c) { return std::tolower(c); } );
783 if (_impl->version == 0) {
784 AliasMap &aliases = *_impl->schema.getAliasMap();
785 for (
auto iter = _impl->asList().begin();
iter != _impl->asList().end(); ++
iter) {
787 if (flagPos != std::string::npos) {
796 ttype.
replace(flagPos, 5,
"flag");
803 }
else if (isInstFlux(*
iter)) {
805 if (endswith(
iter->ttype,
"_err")) {
806 aliases.
set(replace(
iter->ttype,
"_err",
"_instFluxErr"),
iter->ttype);
808 aliases.
set(
iter->ttype +
"_instFlux",
iter->ttype);
810 }
else if (endswith(
iter->ttype,
"_err")) {
816 if (
iter->tccls ==
"Covariance(Point)") {
817 aliases.
set(replace(
iter->ttype,
"_err",
"_yErr"),
iter->ttype +
"_yErr");
818 aliases.
set(replace(
iter->ttype,
"_err",
"_xErr"),
iter->ttype +
"_xErr");
819 aliases.
set(replace(
iter->ttype,
"_err",
"_x_y_Cov"),
iter->ttype +
"_x_y_Cov");
820 }
else if (
iter->tccls ==
"Covariance(Moments)") {
821 aliases.
set(replace(
iter->ttype,
"_err",
"_xxErr"),
iter->ttype +
"_xxErr");
822 aliases.
set(replace(
iter->ttype,
"_err",
"_yyErr"),
iter->ttype +
"_yyErr");
823 aliases.
set(replace(
iter->ttype,
"_err",
"_xyErr"),
iter->ttype +
"_xyErr");
824 aliases.
set(replace(
iter->ttype,
"_err",
"_xx_yy_Cov"),
iter->ttype +
"_xx_yy_Cov");
825 aliases.
set(replace(
iter->ttype,
"_err",
"_xx_xy_Cov"),
iter->ttype +
"_xx_xy_Cov");
826 aliases.
set(replace(
iter->ttype,
"_err",
"_yy_xy_Cov"),
iter->ttype +
"_yy_xy_Cov");
830 }
else if (_impl->version < 3) {
834 AliasMap &aliases = *_impl->schema.getAliasMap();
835 for (
auto iter = _impl->asList().begin();
iter != _impl->asList().end(); ++
iter) {
837 if (_impl->version < 2 && endswith(
name,
"Sigma")) {
840 if (_impl->version < 3 && isInstFlux(*
iter)) {
848 for (
auto iter = _impl->asList().begin();
iter != _impl->asList().end(); ++
iter) {
852 _impl->readers.push_back(
std::move(reader));
854 LOGLS_WARN(
"afw.FitsSchemaInputMapper",
"Format " <<
iter->tform <<
" for column "
856 <<
" not supported; skipping.");
861 (
boost::format(
"Flag field '%s' is is in bit %d (0-indexed) of only %d") %
862 iter->ttype %
iter->bit % _impl->flagKeys.size())
865 _impl->flagKeys[
iter->bit] = _impl->schema.addField<Flag>(
iter->ttype,
iter->doc);
868 _impl->asList().clear();
869 if (_impl->schema.getRecordSize() <= 0) {
872 (
boost::format(
"Non-positive record size: %d; file is corrupt or invalid.") %
873 _impl->schema.getRecordSize()).str()
877 return _impl->schema;
881 if (!_impl->flagKeys.empty()) {
882 fits.readTableArray<
bool>(
row, _impl->flagColumn, _impl->flagKeys.size(), _impl->flagWorkspace.get());
883 for (
std::size_t bit = 0; bit < _impl->flagKeys.size(); ++bit) {
884 record.
set(_impl->flagKeys[bit], _impl->flagWorkspace[bit]);
887 if (_impl->nRowsToPrep != 1 &&
row % _impl->nRowsToPrep == 0) {
891 for (
auto & reader : _impl->readers) {
892 reader->prepRead(
row, size,
fits);
895 for (
auto const & reader : _impl->readers) {
896 reader->readCell(record,
row,
fits, _impl->archive);
table::Key< std::string > name
Key< Flag > const & target
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
LSST DM logging module built on log4cxx.
#define LOGLS_WARN(logger, message)
Log a warn-level message using an iostream-based interface.
An exception thrown when problems are found when reading or writing FITS files.
A simple struct that combines the two arguments that must be passed to most cfitsio routines and cont...
Mapping class that holds aliases for a Schema.
void set(std::string const &alias, std::string const &target)
Add an alias to the schema or replace an existing one.
Base class for all records.
void set(Key< T > const &key, U const &value)
Set value of a field for the given key.
Defines the fields and offsets for a table.
Polymorphic reader interface used to read different kinds of objects from one or more FITS binary tab...
Class for storing ordered metadata with comments.
T get(std::string const &name) const
Get the last value for a property name (possibly hierarchical).
std::string const & getComment(std::string const &name) const
Get the comment for a string property name (possibly hierarchical).
std::vector< T > getArray(std::string const &name) const
Get the vector of values for a property name (possibly hierarchical).
std::vector< std::string > getOrderedNames() const
Get the list of property names, in the order they were added.
virtual void remove(std::string const &name)
Remove all values for a property name (possibly hierarchical).
bool exists(std::string const &name) const
Determine if a name (possibly hierarchical) exists.
A coordinate class intended to represent absolute positions (2-d specialization).
Point in an unspecified spherical coordinate system.
Reports attempts to exceed implementation-defined length limits for some classes.
Reports attempts to access elements using an invalid key.
int computeCovariancePackedSize(int size)
Defines the packed size of a covariance matrices.
int indexCovariance(int i, int j)
Defines the ordering of packed covariance matrices.
CoordinateType
Enum used to set units for geometric FunctorKeys.
constexpr AngleUnit radians
constant with units of radians
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
A base class for image defects.
table::Key< table::Array< int > > _size
Field base class default implementation (used for numeric scalars and lsst::geom::Angle).
A structure that describes a field as a collection of related strings read from the FITS header.