LSSTApplications  18.1.0
LSSTDataManagementBasePackage
Source.cc
Go to the documentation of this file.
1 // -*- lsst-c++ -*-
2 #include <typeinfo>
3 
4 #include "boost/iterator/transform_iterator.hpp"
5 
9 #include "lsst/afw/geom/SkyWcs.h"
14 
15 namespace lsst {
16 namespace afw {
17 namespace table {
18 
19 //-----------------------------------------------------------------------------------------------------------
20 //----- PersistenceHelpers ----------------------------------------------------------------------------------
21 //-----------------------------------------------------------------------------------------------------------
22 
23 namespace {
24 
25 struct PersistenceHelper {
26  SchemaMapper mapper;
27  afw::table::Key<int> footprintKey;
28 };
29 
30 } // namespace
31 
32 //-----------------------------------------------------------------------------------------------------------
33 //----- SourceFitsWriter ------------------------------------------------------------------------------------
34 //-----------------------------------------------------------------------------------------------------------
35 
36 // A custom FitsWriter for Sources. It also sets the AFW_TYPE key to SOURCE, which should
37 // ensure we use SourceFitsReader to read it.
38 
39 // Because it also holds Footprints, a SourceCatalog isn't persisted with the same Schema
40 // it has in-memory; instead, it's saved in a Schema that has an additional int field
41 // appended to the end, which contains the afw::table::io "archive ID" that's used to
42 // extract the Footprint from additional FITS HDUs. (If we disable saving Footprints via
43 // SourceFitsFlags, we do save using the original Schema).
44 
45 // The only public access point to this class is SourceTable::makeFitsWriter. If we
46 // subclass SourceTable someday, it may be necessary to put SourceFitsWriter in a header
47 // file so we can subclass it too.
48 
49 namespace {
50 
51 class SourceFitsWriter : public io::FitsWriter {
52 public:
53  explicit SourceFitsWriter(Fits *fits, int flags) : io::FitsWriter(fits, flags) {}
54 
55 protected:
56  void _writeTable(std::shared_ptr<BaseTable const> const &table, std::size_t nRows) override;
57 
58  void _writeRecord(BaseRecord const &record) override;
59 
60  void _finish() override {
61  if (!(_flags & SOURCE_IO_NO_FOOTPRINTS)) {
62  _archive.writeFits(*_fits);
63  }
64  }
65 
66 private:
67  SchemaMapper _mapper;
68  std::shared_ptr<BaseRecord> _outRecord;
70  Key<int> _footprintKey;
71  io::OutputArchive _archive;
72 };
73 
74 void SourceFitsWriter::_writeTable(std::shared_ptr<BaseTable const> const &t, std::size_t nRows) {
76  if (!table) {
78  "Cannot use a SourceFitsWriter on a non-Source table.");
79  }
80  if (!(_flags & SOURCE_IO_NO_FOOTPRINTS)) {
81  _mapper = SchemaMapper(t->getSchema(), true);
82  _mapper.addMinimalSchema(t->getSchema(), true);
83  _footprintKey = _mapper.editOutputSchema().addField<int>("footprint", "archive ID for Footprint");
84  _outTable = BaseTable::make(_mapper.getOutputSchema());
85  std::shared_ptr<daf::base::PropertyList> metadata = table->getMetadata();
86  if (metadata) {
87  metadata = std::static_pointer_cast<daf::base::PropertyList>(metadata->deepCopy());
88  } else {
89  metadata.reset(new daf::base::PropertyList());
90  }
91  // HDU 0 is empty (primary HDU can't be a table)
92  // HDU 1 is the SourceCatalog's records
93  // HDU 2 is the index for the afw::table::io archive that holds more complex objects
94  //
95  // Historically the AR_HDU keyword was 1-indexed (see RFC-304), and to maintain file compatibility
96  // this is still the case so we're setting AR_HDU to 3 == 2 + 1
97  metadata->set("AR_HDU", 3,
98  "HDU (1-indexed) containing the archive index for non-record data (e.g. Footprints)");
99  _outTable->setMetadata(metadata);
100  _outRecord = _outTable->makeRecord(); // make temporary record to use as a workspace
101  io::FitsWriter::_writeTable(_outTable, nRows);
102  } else {
103  io::FitsWriter::_writeTable(table, nRows);
104  }
105  _fits->writeKey("AFW_TYPE", "SOURCE", "Tells lsst::afw to load this as a Source table.");
106 }
107 
108 void SourceFitsWriter::_writeRecord(BaseRecord const &r) {
109  SourceRecord const &record = static_cast<SourceRecord const &>(r);
110  if (!(_flags & SOURCE_IO_NO_FOOTPRINTS)) {
111  _outRecord->assign(record, _mapper);
113  if (footprint) {
114  if ((_flags & SOURCE_IO_NO_HEAVY_FOOTPRINTS) && footprint->isHeavy()) {
115  footprint.reset(new afw::detection::Footprint(*footprint));
116  }
117  int footprintArchiveId = _archive.put(footprint);
118  _outRecord->set(_footprintKey, footprintArchiveId);
119  }
120  io::FitsWriter::_writeRecord(*_outRecord);
121  } else {
123  }
124 }
125 
126 } // namespace
127 
128 //-----------------------------------------------------------------------------------------------------------
129 //----- SourceFitsReader ------------------------------------------------------------------------------------
130 //-----------------------------------------------------------------------------------------------------------
131 
132 // A custom FitsReader for Sources - this reads footprints in addition to the regular fields. It
133 // gets registered with name SOURCE, so it should get used whenever we read a table with AFW_TYPE
134 // set to that value. (The actual SourceFitsReader class is a bit further down, after some helper
135 // classes.)
136 
137 // As noted in the comments for SourceFitsWriter, we add a column for the Footprint archive ID when
138 // we save a SourceCatalog.
139 
140 // Things are a bit more complicated than that when reading, because we also need to be able to read files
141 // saved with an older version of the pipeline, in which there were 2-5 additional columns, all variable-
142 // length arrays, holding the Spans, Peaks, and HeavyFootprint arrays. Those are handled by explicit
143 // calls to the FITS I/O routines here.
144 
145 // The only public access point to this class is through the registry. If we subclass SourceTable
146 // someday, it may be necessary to put SourceFitsReader in a header file so we can subclass it too.
147 
148 namespace {
149 
150 // FitsColumnReader subclass for backwards-compatible Footprint reading from variable-length arrays
151 class OldSourceFootprintReader : public io::FitsColumnReader {
152 public:
153  static int readSpecialColumn(io::FitsSchemaInputMapper &mapper, daf::base::PropertyList &metadata,
154  bool stripMetadata, std::string const &name) {
155  int column = metadata.get(name, 0);
156  --column; // switch from 1-indexed to 0-indexed convention
157  if (column >= 0) {
158  if (stripMetadata) {
159  metadata.remove(name);
160  }
161  mapper.erase(name);
162  }
163  return column;
164  }
165 
166  static void setup(io::FitsSchemaInputMapper &mapper, daf::base::PropertyList &metadata, int ioFlags,
167  bool stripMetadata) {
168  std::unique_ptr<OldSourceFootprintReader> reader(new OldSourceFootprintReader());
169  reader->_spanCol = readSpecialColumn(mapper, metadata, stripMetadata, "SPANCOL");
170  reader->_peakCol = readSpecialColumn(mapper, metadata, stripMetadata, "PEAKCOL");
171  reader->_heavyPixCol = readSpecialColumn(mapper, metadata, stripMetadata, "HVYPIXCO");
172  reader->_heavyMaskCol = readSpecialColumn(mapper, metadata, stripMetadata, "HVYMSKCO");
173  reader->_heavyVarCol = readSpecialColumn(mapper, metadata, stripMetadata, "HVYVARCO");
174  if ((ioFlags & SOURCE_IO_NO_FOOTPRINTS) || mapper.hasArchive()) {
175  return; // don't want to load anything, so we're done after just removing the special columns
176  }
177  if (ioFlags & SOURCE_IO_NO_HEAVY_FOOTPRINTS) {
178  reader->_heavyPixCol = -1;
179  reader->_heavyMaskCol = -1;
180  reader->_heavyVarCol = -1;
181  }
182  // These checks are really basically assertions - they should only happen if we get
183  // a corrupted catalog - but we still don't want to crash if that happens.
184  if ((reader->_spanCol >= 0) != (reader->_peakCol >= 0)) {
185  throw LSST_EXCEPT(afw::fits::FitsError,
186  "Corrupted catalog: either both or none of the Footprint Span/Peak columns "
187  "must be present.");
188  }
189  if (reader->_spanCol < 0) {
190  return;
191  }
192  if ((reader->_heavyPixCol >= 0) != (reader->_heavyMaskCol >= 0) ||
193  (reader->_heavyPixCol >= 0) != (reader->_heavyVarCol >= 0)) {
194  throw LSST_EXCEPT(
195  afw::fits::FitsError,
196  "Corrupted catalog: either all or none of the HeavyFootprint columns must be present.");
197  }
198  if (reader->_heavyPixCol >= 0 && reader->_spanCol < 0) {
199  throw LSST_EXCEPT(afw::fits::FitsError,
200  "Corrupted catalog: HeavyFootprint columns with no Span/Peak columns.");
201  }
202  // If we do want to load old-style Footprints, add the column reader to the mapper.
203  mapper.customize(std::move(reader));
204  }
205 
206  void readCell(BaseRecord &baseRecord, std::size_t row, fits::Fits &fits,
207  std::shared_ptr<io::InputArchive> const &archive) const override {
208  SourceRecord &record = static_cast<SourceRecord &>(baseRecord);
209  std::vector<geom::Span> spansVector;
210 
211  // Load a regular Footprint from the span and peak columns.
212  int spanElementCount = fits.getTableArraySize(row, _spanCol);
213  int peakElementCount = fits.getTableArraySize(row, _peakCol);
214  if (spanElementCount) {
215  if (spanElementCount % 3) {
216  throw LSST_EXCEPT(
217  afw::fits::FitsError,
219  fits.fptr, fits.status,
220  boost::format("Number of span elements (%d) must divisible by 3 (row %d)") %
221  spanElementCount % row));
222  }
223  std::vector<int> spanElements(spanElementCount);
224  fits.readTableArray(row, _spanCol, spanElementCount, &spanElements.front());
225  std::vector<int>::iterator j = spanElements.begin();
226  while (j != spanElements.end()) {
227  int y = *j++;
228  int x0 = *j++;
229  int x1 = *j++;
230  spansVector.push_back(geom::Span(y, x0, x1));
231  }
232  }
233  std::shared_ptr<Footprint> fp = std::make_shared<detection::Footprint>(
234  std::make_shared<geom::SpanSet>(std::move(spansVector)));
235  if (peakElementCount) {
236  if (peakElementCount % 3) {
237  throw LSST_EXCEPT(
238  afw::fits::FitsError,
240  fits.fptr, fits.status,
241  boost::format("Number of peak elements (%d) must divisible by 3 (row %d)") %
242  peakElementCount % row));
243  }
244  std::vector<float> peakElements(peakElementCount);
245  fits.readTableArray(row, _peakCol, peakElementCount, &peakElements.front());
246  std::vector<float>::iterator j = peakElements.begin();
247  while (j != peakElements.end()) {
248  float x = *j++;
249  float y = *j++;
250  float value = *j++;
251  fp->addPeak(x, y, value);
252  }
253  }
254  record.setFootprint(fp);
255 
256  // If we're setup to read HeavyFootprints
257  if (_heavyPixCol < 0) {
258  return;
259  }
260  int heavyPixElementCount = fits.getTableArraySize(row, _heavyPixCol);
261  int heavyMaskElementCount = fits.getTableArraySize(row, _heavyMaskCol);
262  int heavyVarElementCount = fits.getTableArraySize(row, _heavyVarCol);
263  if (heavyPixElementCount > 0) {
264  int N = fp->getArea();
265  if ((heavyPixElementCount != N) || (heavyMaskElementCount != N) || (heavyVarElementCount != N)) {
266  throw LSST_EXCEPT(
267  afw::fits::FitsError,
269  fits.fptr, fits.status,
270  boost::format("Number of HeavyFootprint elements (pix %d, mask %d, var %d) "
271  "must all be equal to footprint area (%d)") %
272  heavyPixElementCount % heavyMaskElementCount % heavyVarElementCount %
273  N));
274  }
275  // float HeavyFootprints were the only kind we ever saved using the old format
276  typedef detection::HeavyFootprint<float, image::MaskPixel, image::VariancePixel> HeavyFootprint;
277  std::shared_ptr<HeavyFootprint> heavy = std::make_shared<HeavyFootprint>(*fp);
278  fits.readTableArray(row, _heavyPixCol, N, heavy->getImageArray().getData());
279  fits.readTableArray(row, _heavyMaskCol, N, heavy->getMaskArray().getData());
280  fits.readTableArray(row, _heavyVarCol, N, heavy->getVarianceArray().getData());
281  record.setFootprint(heavy);
282  }
283  }
284 
285 private:
286  int _spanCol;
287  int _peakCol;
288  int _heavyPixCol;
289  int _heavyMaskCol;
290  int _heavyVarCol;
291 };
292 
293 // FitsColumnReader for new-style Footprint persistence using archives.
294 class SourceFootprintReader : public io::FitsColumnReader {
295 public:
296  static void setup(io::FitsSchemaInputMapper &mapper, int ioFlags) {
297  auto item = mapper.find("footprint");
298  if (item) {
299  if (mapper.hasArchive()) {
301  new SourceFootprintReader(ioFlags & SOURCE_IO_NO_HEAVY_FOOTPRINTS, item->column));
302  mapper.customize(std::move(reader));
303  }
304  mapper.erase(item);
305  }
306  }
307 
308  SourceFootprintReader(bool noHeavy, int column) : _noHeavy(noHeavy), _column(column) {}
309 
310  void readCell(BaseRecord &record, std::size_t row, fits::Fits &fits,
311  std::shared_ptr<io::InputArchive> const &archive) const override {
312  int id = 0;
313  fits.readTableScalar<int>(row, _column, id);
315  if (_noHeavy && footprint->isHeavy()) {
316  // It sort of defeats the purpose of the flag if we have to do the I/O to read
317  // a HeavyFootprint before we can downgrade it to a regular Footprint, but that's
318  // what we're going to do - at least this will save on on some memory usage, which
319  // might still be useful. It'd be really hard to fix this
320  // (because we have no way to pass something like the ioFlags to the InputArchive).
321  // The good news is that if someone's concerned about performance of reading
322  // SourceCatalogs, they'll almost certainly use SOURCE_IO_NO_FOOTPRINTS, which
323  // will do what we want. SOURCE_IO_NO_HEAVY_FOOTPRINTS is more useful for writing
324  // sources, and that still works just fine.
325  footprint.reset(new Footprint(*footprint));
326  }
327  static_cast<SourceRecord &>(record).setFootprint(footprint);
328  }
329 
330 private:
331  bool _noHeavy;
332  int _column;
333 };
334 
335 class SourceFitsReader : public io::FitsReader {
336 public:
337  SourceFitsReader() : afw::table::io::FitsReader("SOURCE") {}
338 
339  std::shared_ptr<BaseTable> makeTable(io::FitsSchemaInputMapper &mapper,
340  std::shared_ptr<daf::base::PropertyList> metadata, int ioFlags,
341  bool stripMetadata) const override {
342  // Look for old-style persistence of Footprints. If we have both that and an archive, we
343  // load the footprints from the archive, but still need to remove the old-style header keys
344  // from the metadata and the corresponding fields from the FitsSchemaInputMapper.
345  OldSourceFootprintReader::setup(mapper, *metadata, ioFlags, stripMetadata);
346  // Look for new-style persistence of Footprints. We'll only read them if we have an archive,
347  // but we'll strip fields out regardless.
348  SourceFootprintReader::setup(mapper, ioFlags);
349  std::shared_ptr<SourceTable> table = SourceTable::make(mapper.finalize());
350  table->setMetadata(metadata);
351  return table;
352  }
353 
354  bool usesArchive(int ioFlags) const override { return !(ioFlags & SOURCE_IO_NO_FOOTPRINTS); }
355 };
356 
357 // registers the reader so FitsReader::make can use it.
358 static SourceFitsReader const sourceFitsReader;
359 
360 } // namespace
361 
362 //-----------------------------------------------------------------------------------------------------------
363 //----- SourceTable/Record member function implementations --------------------------------------------------
364 //-----------------------------------------------------------------------------------------------------------
365 
366 SourceRecord::~SourceRecord() = default;
367 
368 void SourceRecord::updateCoord(geom::SkyWcs const &wcs) { setCoord(wcs.pixelToSky(getCentroid())); }
369 
371  setCoord(wcs.pixelToSky(get(key)));
372 }
373 
375  try {
376  SourceRecord const &s = dynamic_cast<SourceRecord const &>(other);
377  _footprint = s._footprint;
378  } catch (std::bad_cast &) {
379  }
380 }
381 
384  if (!checkSchema(schema)) {
386  "Schema for Source must contain at least the keys defined by getMinimalSchema().");
387  }
388  return std::shared_ptr<SourceTable>(new SourceTable(schema, idFactory));
389 }
390 
392  : SimpleTable(schema, idFactory), _slots(schema) {}
393 
394 SourceTable::SourceTable(SourceTable const &other) : SimpleTable(other), _slots(other._slots) {}
395 // Delegate to copy constructor for backward compatibility
397 
399  if (alias.compare(0, 4, "slot") != 0) {
400  return;
401  }
402  _slots.handleAliasChange(alias, getSchema());
403 }
404 
405 SourceTable::MinimalSchema::MinimalSchema() {
407  parent = schema.addField<RecordId>("parent", "unique ID of parent source");
408  schema.getCitizen().markPersistent();
409 }
410 
411 SourceTable::MinimalSchema &SourceTable::getMinimalSchema() {
412  static MinimalSchema it;
413  return it;
414 }
415 
416 std::shared_ptr<io::FitsWriter> SourceTable::makeFitsWriter(fits::Fits *fitsfile, int flags) const {
417  return std::make_shared<SourceFitsWriter>(fitsfile, flags);
418 }
419 
421  return std::shared_ptr<SourceTable>(new SourceTable(*this));
422 }
423 
425  auto record = constructRecord<SourceRecord>();
426  if (getIdFactory()) record->setId((*getIdFactory())());
427  return record;
428 }
429 
430 template class CatalogT<SourceRecord>;
431 template class CatalogT<SourceRecord const>;
432 
433 template class SortedCatalogT<SourceRecord>;
435 } // namespace table
436 } // namespace afw
437 } // namespace lsst
void handleAliasChange(std::string const &alias) override
Definition: Source.cc:398
Defines the fields and offsets for a table.
Definition: Schema.h:50
A 2-dimensional celestial WCS that transform pixels to ICRS RA/Dec, using the LSST standard for pixel...
Definition: SkyWcs.h:117
static Schema makeMinimalSchema()
Return a minimal schema for Simple tables and records.
Definition: Simple.h:140
void updateCoord(geom::SkyWcs const &wcs)
Update the coord field using the given Wcs and the field in the centroid slot.
Definition: Source.cc:368
T front(T... args)
A custom container class for records, based on std::vector.
Definition: Catalog.h:97
lsst::geom::SpherePoint pixelToSky(lsst::geom::Point2D const &pixel) const
Compute sky position(s) from pixel position(s)
Definition: SkyWcs.h:334
Schema getSchema() const
Return the table&#39;s schema.
Definition: BaseTable.h:137
int y
Definition: SpanSet.cc:49
T end(T... args)
table::Key< int > id
Definition: Detector.cc:166
STL class.
A simple struct that combines the two arguments that must be passed to most cfitsio routines and cont...
Definition: fits.h:297
virtual void _assign(BaseRecord const &other)
Called by assign() after transferring fields to allow subclass data members to be copied...
Definition: Source.cc:374
table::Key< table::Array< std::uint8_t > > wcs
Definition: SkyWcs.cc:71
Fits * fits
Definition: FitsWriter.cc:90
STL class.
std::shared_ptr< BaseRecord > _makeRecord() override
Default-construct an associated record (protected implementation).
Definition: Source.cc:424
T push_back(T... args)
Table class that must contain a unique ID field and a celestial coordinate field. ...
Definition: Simple.h:102
A base class for image defects.
Custom catalog class for record/table subclasses that are guaranteed to have an ID, and should generally be sorted by that ID.
Definition: fwd.h:63
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:168
table::Schema schema
Definition: Camera.cc:161
T reset(T... args)
SchemaMapper _mapper
Definition: Exposure.cc:220
T dynamic_pointer_cast(T... args)
Reports errors in the logical structure of the program.
Definition: Runtime.h:46
solver_t * s
T move(T... args)
virtual void _writeTable(std::shared_ptr< BaseTable const > const &table, std::size_t nRows)
Write a table and its schema.
Definition: FitsWriter.cc:103
double x
std::shared_ptr< BaseTable > _clone() const override
Clone implementation with noncovariant return types.
Definition: Source.cc:420
std::shared_ptr< io::OutputArchive > _archive
Definition: Exposure.cc:217
static std::shared_ptr< SourceTable > make(Schema const &schema, std::shared_ptr< IdFactory > const &idFactory)
Construct a new table.
Definition: Source.cc:382
T get(T... args)
Table class that contains measurements made on a single exposure.
Definition: Source.h:221
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
STL class.
STL class.
Base class for all records.
Definition: BaseRecord.h:31
afw::table::Key< int > footprintKey
Definition: Source.cc:27
static std::shared_ptr< BaseTable > make(Schema const &schema)
Construct a new table.
Definition: BaseTable.cc:121
T begin(T... args)
Key< U > key
Definition: Schema.cc:281
Read/write heavy footprints as non-heavy footprints.
Definition: Source.h:58
lsst::afw::detection::Footprint Footprint
Definition: Source.h:61
SchemaMapper mapper
Definition: Source.cc:26
virtual void _writeRecord(BaseRecord const &source)
Write an individual record.
Definition: FitsWriter.cc:189
Reports invalid arguments.
Definition: Runtime.h:66
ItemVariant const * other
Definition: Schema.cc:56
Do not read/write footprints at all.
Definition: Source.h:57
Record class that contains measurements made on a single exposure.
Definition: Source.h:82
std::shared_ptr< IdFactory > getIdFactory()
Return the object that generates IDs for the table (may be null).
Definition: Simple.h:155
std::string makeErrorMessage(std::string const &fileName="", int status=0, std::string const &msg="")
Return an error message reflecting FITS I/O errors.
Definition: fits.cc:407
SourceTable(Schema const &schema, std::shared_ptr< IdFactory > const &idFactory)
Definition: Source.cc:391
T compare(T... args)
int row
Definition: CR.cc:145
void handleAliasChange(std::string const &alias, Schema const &schema)
Handle a callback from an AliasMap informing the table that an alias has changed. ...
Definition: slots.cc:116