LSST Applications g11f7dcd041+3309dba821,g1cd03abf6b+5538b8be6d,g1ce3e0751c+f991eae79d,g28da252d5a+9cd75b2bbf,g2bbee38e9b+b474e213a6,g2bc492864f+b474e213a6,g2cdde0e794+1e3907a55c,g2e95ee2dae+b3f7f8b914,g347aa1857d+b474e213a6,g35bb328faa+b86e4b8053,g3a166c0a6a+b474e213a6,g4322eb9e3a+f7440700bb,g461a3dce89+b86e4b8053,g52b1c1532d+b86e4b8053,g78056777b3+b52dbcd3a1,g858d7b2824+d6777fb84b,g8cd86fa7b1+15b40b54d8,g9ddcbc5298+f24b38b85a,ga61afda322+d6777fb84b,ga68bd05109+bf27b4c466,gae0086650b+b86e4b8053,gbb886bcc26+43b45a029f,gbd462c55f0+538be87034,gc28159a63d+b474e213a6,gc2a6998b7e+f95f64aeae,gc30aee3386+6076c0b52b,gcaf7e4fdec+d6777fb84b,gcd45df26be+d6777fb84b,gcdd4ae20e8+490f2fb51f,gce08ada175+bde66f6e3f,gcf0d15dbbd+490f2fb51f,gdaeeff99f8+006e14e809,gdbce86181e+61a8f2cffc,ge3d4d395c2+5a6209b744,ge77125146a+06b6d19373,ge79ae78c31+b474e213a6,gf048a9a2f4+060b48ea8f,gf9c3535cb8+b86e4b8053,w.2024.28
LSST Data Management Base Package
Loading...
Searching...
No Matches
Source.cc
Go to the documentation of this file.
1// -*- lsst-c++ -*-
2#include <typeinfo>
3
4#include "boost/format.hpp"
5
14
15namespace lsst {
16namespace afw {
17namespace table {
18
19//-----------------------------------------------------------------------------------------------------------
20//----- PersistenceHelpers ----------------------------------------------------------------------------------
21//-----------------------------------------------------------------------------------------------------------
22
23namespace {
24
25struct 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
49namespace {
50
51class SourceFitsWriter : public io::FitsWriter {
52public:
53 explicit SourceFitsWriter(Fits *fits, int flags) : io::FitsWriter(fits, flags) {}
54
55protected:
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
66private:
67 SchemaMapper _mapper;
70 Key<int> _footprintKey;
71 io::OutputArchive _archive;
72};
73
74void SourceFitsWriter::_writeTable(std::shared_ptr<BaseTable const> const &t, std::size_t nRows) {
75 std::shared_ptr<SourceTable const> table = std::dynamic_pointer_cast<SourceTable const>(t);
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
108void 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);
112 std::shared_ptr<afw::detection::Footprint> footprint = record.getFootprint();
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 {
122 io::FitsWriter::_writeRecord(record);
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
148namespace {
149
150// FitsColumnReader subclass for backwards-compatible Footprint reading from variable-length arrays
151class OldSourceFootprintReader : public io::FitsColumnReader {
152public:
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 std::size_t spanElementCount = fits.getTableArraySize(row, _spanCol);
213 std::size_t peakElementCount = fits.getTableArraySize(row, _peakCol);
214 if (spanElementCount) {
215 if (spanElementCount % 3) {
216 throw LSST_EXCEPT(
217 afw::fits::FitsError,
218 afw::fits::makeErrorMessage(
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.emplace_back(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,
239 afw::fits::makeErrorMessage(
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 std::size_t heavyPixElementCount = fits.getTableArraySize(row, _heavyPixCol);
261 std::size_t heavyMaskElementCount = fits.getTableArraySize(row, _heavyMaskCol);
262 std::size_t heavyVarElementCount = fits.getTableArraySize(row, _heavyVarCol);
263 if (heavyPixElementCount > 0) {
264 std::size_t N = fp->getArea();
265 if ((heavyPixElementCount != N) || (heavyMaskElementCount != N) || (heavyVarElementCount != N)) {
266 throw LSST_EXCEPT(
267 afw::fits::FitsError,
268 afw::fits::makeErrorMessage(
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 using HeavyFootprint = detection::HeavyFootprint<float>;
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
285private:
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.
294class SourceFootprintReader : public io::FitsColumnReader {
295public:
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);
314 std::shared_ptr<Footprint> footprint = archive->get<Footprint>(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
330private:
331 bool _noHeavy;
332 int _column;
333};
334
335class SourceFitsReader : public io::FitsReader {
336public:
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.
358static SourceFitsReader const sourceFitsReader;
359
360} // namespace
361
362//-----------------------------------------------------------------------------------------------------------
363//----- SourceTable/Record member function implementations --------------------------------------------------
364//-----------------------------------------------------------------------------------------------------------
365
367
368void SourceRecord::updateCoord(geom::SkyWcs const &wcs, bool include_covariance) {
369 lsst::geom::Point2D center = getCentroid();
370 setCoord(wcs.pixelToSky(center));
371 if (include_covariance) {
372 // Get coordinate covariance:
373 auto err = getCentroidErr();
375 Eigen::Matrix2f skyCov = calculateCoordCovariance(wcs, center, err);
377 }
378}
379
380void SourceRecord::updateCoord(geom::SkyWcs const &wcs, PointKey<double> const &key, bool include_covariance) {
381 lsst::geom::Point2D center = get(key);
382 setCoord(wcs.pixelToSky(center));
383 if (include_covariance) {
384 // Get coordinate covariance:
385 auto err = getCentroidErr();
387 Eigen::Matrix2f skyCov = calculateCoordCovariance(wcs, center, err);
389 }
390}
391
393 try {
394 SourceRecord const &s = dynamic_cast<SourceRecord const &>(other);
395 _footprint = s._footprint;
396 } catch (std::bad_cast &) {
397 }
398}
399
401 std::shared_ptr<IdFactory> const &idFactory) {
402 if (!checkSchema(schema)) {
404 "Schema for Source must contain at least the keys defined by getMinimalSchema().");
405 }
406 return std::shared_ptr<SourceTable>(new SourceTable(schema, idFactory));
407}
408
410 : SimpleTable(schema, idFactory), _slots(schema) {}
411
412SourceTable::SourceTable(SourceTable const &other) = default;
413// Delegate to copy constructor for backward compatibility
415
417 if (alias.compare(0, 4, "slot") != 0) {
418 return;
419 }
421}
422
423SourceTable::MinimalSchema::MinimalSchema() {
425 parent = schema.addField<RecordId>("parent", "unique ID of parent source");
426}
427
428SourceTable::MinimalSchema &SourceTable::getMinimalSchema() {
429 static MinimalSchema it;
430 return it;
431}
432
434 return std::make_shared<SourceFitsWriter>(fitsfile, flags);
435}
436
440
442 auto record = constructRecord<SourceRecord>();
443 if (getIdFactory()) record->setId((*getIdFactory())());
444 return record;
445}
446
447template class CatalogT<SourceRecord>;
448template class CatalogT<SourceRecord const>;
449
450template class SortedCatalogT<SourceRecord>;
452} // namespace table
453} // namespace afw
454} // namespace lsst
table::Key< int > id
Definition Detector.cc:162
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
SchemaMapper * mapper
int y
Definition SpanSet.cc:48
afw::table::Key< int > footprintKey
Definition Source.cc:27
A simple struct that combines the two arguments that must be passed to most cfitsio routines and cont...
Definition fits.h:308
Tag types used to declare specialized field types.
Definition misc.h:31
Base class for all records.
Definition BaseRecord.h:31
Schema getSchema() const
Return the table's schema.
Definition BaseTable.h:137
static ErrorKey getErrorKey(Schema const &schema)
Defines the fields and offsets for a table.
Definition Schema.h:51
Table class that must contain a unique ID field and a celestial coordinate field.
Definition Simple.h:102
static Schema makeMinimalSchema()
Return a minimal schema for Simple tables and records.
Definition Simple.h:140
std::shared_ptr< IdFactory > getIdFactory()
Return the object that generates IDs for the table (may be null).
Definition Simple.h:155
Record class that contains measurements made on a single exposure.
Definition Source.h:78
virtual void _assign(BaseRecord const &other)
Called by assign() after transferring fields to allow subclass data members to be copied.
Definition Source.cc:392
void updateCoord(geom::SkyWcs const &wcs, bool include_covariance=true)
Update the coord field using the given Wcs and the field in the centroid slot.
Definition Source.cc:368
Table class that contains measurements made on a single exposure.
Definition Source.h:217
static std::shared_ptr< SourceTable > make(Schema const &schema, std::shared_ptr< IdFactory > const &idFactory)
Construct a new table.
Definition Source.cc:400
SourceTable(Schema const &schema, std::shared_ptr< IdFactory > const &idFactory)
Definition Source.cc:409
std::shared_ptr< BaseRecord > _makeRecord() override
Default-construct an associated record (protected implementation).
Definition Source.cc:441
std::shared_ptr< io::FitsWriter > makeFitsWriter(fits::Fits *fitsfile, int flags) const override
Definition Source.cc:433
std::shared_ptr< BaseTable > _clone() const override
Clone implementation with noncovariant return types.
Definition Source.cc:437
void handleAliasChange(std::string const &alias) override
Definition Source.cc:416
Reports invalid arguments.
Definition Runtime.h:66
Reports errors in the logical structure of the program.
Definition Runtime.h:46
T emplace_back(T... args)
daf::base::PropertySet * set
Definition fits.cc:931
T get(T... args)
T move(T... args)
Eigen::Matrix2f calculateCoordCovariance(geom::SkyWcs const &wcs, lsst::geom::Point2D center, Eigen::Matrix2f err)
Calculate covariance for sky coordinates.
Definition wcsUtils.cc:95
@ SOURCE_IO_NO_FOOTPRINTS
Do not read/write footprints at all.
Definition Source.h:53
T reset(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
std::shared_ptr< io::OutputArchive > _archive
Definition Exposure.cc:217
SchemaMapper _mapper
Definition Exposure.cc:220