LSST Applications g0b6bd0c080+a72a5dd7e6,g1182afd7b4+2a019aa3bb,g17e5ecfddb+2b8207f7de,g1d67935e3f+06cf436103,g38293774b4+ac198e9f13,g396055baef+6a2097e274,g3b44f30a73+6611e0205b,g480783c3b1+98f8679e14,g48ccf36440+89c08d0516,g4b93dc025c+98f8679e14,g5c4744a4d9+a302e8c7f0,g613e996a0d+e1c447f2e0,g6c8d09e9e7+25247a063c,g7271f0639c+98f8679e14,g7a9cd813b8+124095ede6,g9d27549199+a302e8c7f0,ga1cf026fa3+ac198e9f13,ga32aa97882+7403ac30ac,ga786bb30fb+7a139211af,gaa63f70f4e+9994eb9896,gabf319e997+ade567573c,gba47b54d5d+94dc90c3ea,gbec6a3398f+06cf436103,gc6308e37c7+07dd123edb,gc655b1545f+ade567573c,gcc9029db3c+ab229f5caf,gd01420fc67+06cf436103,gd877ba84e5+06cf436103,gdb4cecd868+6f279b5b48,ge2d134c3d5+cc4dbb2e3f,ge448b5faa6+86d1ceac1d,gecc7e12556+98f8679e14,gf3ee170dca+25247a063c,gf4ac96e456+ade567573c,gf9f5ea5b4d+ac198e9f13,gff490e6085+8c2580be5c,w.2022.27
LSST Data Management Base Package
CoaddBoundedField.cc
Go to the documentation of this file.
1// -*- LSST-C++ -*-
2/*
3 * LSST Data Management System
4 * Copyright 2008-2014, 2010 LSST Corporation.
5 *
6 * This product includes software developed by the
7 * LSST Project (http://www.lsst.org/).
8 *
9 * This program is free software: you can redistribute it and/or modify
10 * it under the terms of the GNU General Public License as published by
11 * the Free Software Foundation, either version 3 of the License, or
12 * (at your option) any later version.
13 *
14 * This program is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 * GNU General Public License for more details.
18 *
19 * You should have received a copy of the LSST License Statement and
20 * the GNU General Public License along with this program. If not,
21 * see <http://www.lsstcorp.org/LegalNotices/>.
22 */
23
25#include "lsst/geom/Box.h"
32
33namespace lsst {
34namespace afw {
35namespace table {
36namespace io {
37
40
41} // namespace io
42} // namespace table
43} // namespace afw
44namespace meas {
45namespace algorithms {
46namespace {
47
48/*
49 * Compare two pointers of the same type
50 *
51 * If both pointers are set then return *a == *b
52 * else return a == b
53 */
54template <typename T>
55bool ptrEquals(T a, T b) {
56 if (a == b) {
57 // test this first for efficiency
58 return true;
59 } else if (a && b) {
60 // both pointers are set, so it is safe to check value equality
61 return *a == *b;
62 }
63 return false;
64}
65
66} // namespace
67
68bool CoaddBoundedFieldElement::operator==(CoaddBoundedFieldElement const& rhs) const {
69 return ptrEquals(field, rhs.field) && ptrEquals(wcs, rhs.wcs) &&
70 ptrEquals(validPolygon, rhs.validPolygon) && (weight == rhs.weight);
71}
72
73CoaddBoundedField::CoaddBoundedField(geom::Box2I const& bbox, std::shared_ptr<afw::geom::SkyWcs const> coaddWcs,
74 ElementVector const& elements)
75 : afw::math::BoundedField(bbox),
76 _throwOnMissing(true),
77 _default(0.0), // unused
78 _coaddWcs(coaddWcs),
79 _elements(elements) {}
80
82 ElementVector const& elements, double default_)
83 : afw::math::BoundedField(bbox),
84 _throwOnMissing(false),
85 _default(default_),
86 _coaddWcs(coaddWcs),
87 _elements(elements) {}
88
89double CoaddBoundedField::evaluate(geom::Point2D const& position) const {
90 auto coord = _coaddWcs->pixelToSky(position);
91 double sum = 0.0;
92 double wSum = 0.0;
93 for (ElementVector::const_iterator i = _elements.begin(); i != _elements.end(); ++i) {
94 geom::Point2D transformedPosition = i->wcs->skyToPixel(coord);
95 bool inValidArea = i->validPolygon ? i->validPolygon->contains(transformedPosition) : true;
96 if (geom::Box2D(i->field->getBBox()).contains(transformedPosition) && inValidArea) {
97 sum += i->weight * i->field->evaluate(transformedPosition);
98 wSum += i->weight;
99 }
100 }
101 if (wSum == 0.0) {
102 if (_throwOnMissing) {
104 (boost::format("No constituent fields to evaluate at point %f, %f") %
105 position.getX() % position.getY())
106 .str());
107 } else {
108 return _default;
109 }
110 }
111 return sum / wSum;
112}
113
114// ---------- Persistence -----------------------------------------------------------------------------------
115
116// For persistence of CoaddBoundedField, we have two catalogs: the first has just one record, and contains
117// the archive ID of the coadd WCS and the parameters that control missing data. The other catalog
118// has one record for each element, with fields corresponding to the data members of the Element struct.
119
120namespace {
121
122// Singleton class that manages the first persistence catalog's schema and keys
123class CoaddBoundedFieldPersistenceKeys1 {
124public:
125 afw::table::Schema schema;
126 afw::table::PointKey<int> bboxMin;
127 afw::table::PointKey<int> bboxMax;
128 afw::table::Key<int> coaddWcs;
129 afw::table::Key<afw::table::Flag> throwOnMissing;
130 afw::table::Key<double> default_;
131
132 static CoaddBoundedFieldPersistenceKeys1 const& get() {
133 static CoaddBoundedFieldPersistenceKeys1 const instance;
134 return instance;
135 }
136
137 // No copying
138 CoaddBoundedFieldPersistenceKeys1(const CoaddBoundedFieldPersistenceKeys1&) = delete;
139 CoaddBoundedFieldPersistenceKeys1& operator=(const CoaddBoundedFieldPersistenceKeys1&) = delete;
140
141 // No moving
142 CoaddBoundedFieldPersistenceKeys1(CoaddBoundedFieldPersistenceKeys1&&) = delete;
143 CoaddBoundedFieldPersistenceKeys1& operator=(CoaddBoundedFieldPersistenceKeys1&&) = delete;
144
145private:
146 CoaddBoundedFieldPersistenceKeys1()
147 : schema(),
148 bboxMin(afw::table::PointKey<int>::addFields(schema, "bbox_min",
149 "lower-left corner of bounding box", "pixel")),
150 bboxMax(afw::table::PointKey<int>::addFields(schema, "bbox_max",
151 "upper-right corner of bounding box", "pixel")),
152 coaddWcs(schema.addField<int>("coaddWcs", "archive ID of the coadd's WCS")),
153 throwOnMissing(schema.addField<afw::table::Flag>(
154 "throwOnMissing", "whether to throw an exception on missing data")),
155 default_(schema.addField<double>("default",
156 "default value to use when throwOnMissing is False")) {}
157};
158
159// Singleton class that manages the second persistence catalog's schema and keys
160class CoaddBoundedFieldPersistenceKeys2 {
161public:
162 afw::table::Schema schema;
163 afw::table::Key<int> field;
164 afw::table::Key<int> wcs;
165 afw::table::Key<int> validPolygon;
166 afw::table::Key<double> weight;
167
168 static CoaddBoundedFieldPersistenceKeys2 const& get() {
169 static CoaddBoundedFieldPersistenceKeys2 const instance;
170 return instance;
171 }
172
173 // No copying
174 CoaddBoundedFieldPersistenceKeys2(const CoaddBoundedFieldPersistenceKeys2&) = delete;
175 CoaddBoundedFieldPersistenceKeys2& operator=(const CoaddBoundedFieldPersistenceKeys2&) = delete;
176
177 // No moving
178 CoaddBoundedFieldPersistenceKeys2(CoaddBoundedFieldPersistenceKeys2&&) = delete;
179 CoaddBoundedFieldPersistenceKeys2& operator=(CoaddBoundedFieldPersistenceKeys2&&) = delete;
180
181private:
182 CoaddBoundedFieldPersistenceKeys2()
183 : schema(),
184 field(schema.addField<int>("field", "archive ID of the BoundedField to be coadded")),
185 wcs(schema.addField<int>("wcs", "archive ID of the Wcs associated with this element")),
186 validPolygon(schema.addField<int>("validPolygon",
187 "archive ID of the Polygon associated with this element")),
188 weight(schema.addField<double>("weight", "weight value for this element")) {}
189};
190
191} // namespace
192
194public:
196 read(InputArchive const& archive, CatalogVector const& catalogs) const {
197 CoaddBoundedFieldPersistenceKeys1 const& keys1 = CoaddBoundedFieldPersistenceKeys1::get();
198 CoaddBoundedFieldPersistenceKeys2 const& keys2 = CoaddBoundedFieldPersistenceKeys2::get();
199 LSST_ARCHIVE_ASSERT(catalogs.size() == 2u);
200 LSST_ARCHIVE_ASSERT(catalogs.front().getSchema() == keys1.schema);
201 LSST_ARCHIVE_ASSERT(catalogs.back().getSchema() == keys2.schema);
202 afw::table::BaseRecord const& record1 = catalogs.front().front();
203 ElementVector elements;
204 elements.reserve(catalogs.back().size());
205 for (afw::table::BaseCatalog::const_iterator i = catalogs.back().begin(); i != catalogs.back().end();
206 ++i) {
207 elements.push_back(Element(archive.get<afw::math::BoundedField>(i->get(keys2.field)),
208 archive.get<afw::geom::SkyWcs>(i->get(keys2.wcs)),
209 archive.get<afw::geom::polygon::Polygon>(i->get(keys2.validPolygon)),
210 i->get(keys2.weight)));
211 }
212 return std::make_shared<CoaddBoundedField>(
213 geom::Box2I(record1.get(keys1.bboxMin), record1.get(keys1.bboxMax)),
214 archive.get<afw::geom::SkyWcs>(record1.get(keys1.coaddWcs)), elements,
215 record1.get(keys1.default_));
216 }
217
218 Factory(std::string const& name) : afw::table::io::PersistableFactory(name) {}
219};
220
221namespace {
222
223std::string getCoaddBoundedFieldPersistenceName() { return "CoaddBoundedField"; }
224
225CoaddBoundedField::Factory registration(getCoaddBoundedFieldPersistenceName());
226
227} // namespace
228
229std::string CoaddBoundedField::getPersistenceName() const { return getCoaddBoundedFieldPersistenceName(); }
230
231std::string CoaddBoundedField::getPythonModule() const { return "lsst.meas.algorithms"; }
232
234 CoaddBoundedFieldPersistenceKeys1 const& keys1 = CoaddBoundedFieldPersistenceKeys1::get();
235 CoaddBoundedFieldPersistenceKeys2 const& keys2 = CoaddBoundedFieldPersistenceKeys2::get();
236 afw::table::BaseCatalog cat1 = handle.makeCatalog(keys1.schema);
238 record1->set(keys1.bboxMin, getBBox().getMin());
239 record1->set(keys1.bboxMax, getBBox().getMax());
240 record1->set(keys1.coaddWcs, handle.put(_coaddWcs));
241 record1->set(keys1.default_, _default);
242 handle.saveCatalog(cat1);
243 afw::table::BaseCatalog cat2 = handle.makeCatalog(keys2.schema);
244 for (ElementVector::const_iterator i = _elements.begin(); i != _elements.end(); ++i) {
246 record2->set(keys2.field, handle.put(i->field));
247 record2->set(keys2.wcs, handle.put(i->wcs));
248 record2->set(keys2.validPolygon, handle.put(i->validPolygon));
249 record2->set(keys2.weight, i->weight);
250 }
251 handle.saveCatalog(cat2);
252}
253
255 throw LSST_EXCEPT(pex::exceptions::NotFoundError, "Scaling of CoaddBoundedField is not implemented");
256}
257
258bool CoaddBoundedField::operator==(BoundedField const& rhs) const {
259 auto rhsCasted = dynamic_cast<CoaddBoundedField const*>(&rhs);
260 if (!rhsCasted) return false;
261
262 return (getBBox() == rhsCasted->getBBox()) && (_default == rhsCasted->_default) &&
263 ptrEquals(_coaddWcs, rhsCasted->_coaddWcs) && (_elements == rhsCasted->_elements);
264}
265
266} // namespace algorithms
267} // namespace meas
268} // namespace lsst
table::Key< std::string > name
Definition: Amplifier.cc:116
AmpInfoBoxKey bbox
Definition: Amplifier.cc:117
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
table::Key< int > b
table::Key< int > a
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
Definition: Persistable.h:48
T back(T... args)
T begin(T... args)
A 2-dimensional celestial WCS that transform pixels to ICRS RA/Dec, using the LSST standard for pixel...
Definition: SkyWcs.h:117
Cartesian polygons.
Definition: Polygon.h:59
An abstract base class for 2-d functions defined on an integer bounding boxes.
Definition: BoundedField.h:55
lsst::geom::Box2I getBBox() const
Return the bounding box that defines the region where the field is valid.
Definition: BoundedField.h:112
Base class for all records.
Definition: BaseRecord.h:31
Field< T >::Value get(Key< T > const &key) const
Return the value of a field for the given key.
Definition: BaseRecord.h:151
Iterator class for CatalogT.
Definition: Catalog.h:40
std::shared_ptr< RecordT > addNew()
Create a new record, add it to the end of the catalog, and return a pointer to it.
Definition: Catalog.h:490
A vector of catalogs used by Persistable.
Definition: CatalogVector.h:29
A multi-catalog archive object used to load table::io::Persistable objects.
Definition: InputArchive.h:31
std::shared_ptr< Persistable > get(int id) const
Load the Persistable with the given ID and return it.
An object passed to Persistable::write to allow it to persist itself.
void saveCatalog(BaseCatalog const &catalog)
Save a catalog in the archive.
BaseCatalog makeCatalog(Schema const &schema)
Return a new, empty catalog with the given schema.
int put(Persistable const *obj, bool permissive=false)
Save an object to the archive and return a unique ID that can be used to retrieve it from an InputArc...
static std::shared_ptr< T > dynamicCast(std::shared_ptr< Persistable > const &ptr)
Dynamically cast a shared_ptr.
Definition: Persistable.cc:18
A base class for factory classes used to reconstruct objects from records.
Definition: Persistable.h:228
A floating-point coordinate rectangle geometry.
Definition: Box.h:413
bool contains(Point2D const &point) const noexcept
Return true if the box contains the point.
Definition: Box.cc:322
An integer coordinate rectangle.
Definition: Box.h:55
virtual std::shared_ptr< afw::table::io::Persistable > read(InputArchive const &archive, CatalogVector const &catalogs) const
Construct a new object from the given InputArchive and vector of catalogs.
CoaddBoundedField(geom::Box2I const &bbox, std::shared_ptr< afw::geom::SkyWcs const > coaddWcs, ElementVector const &elements)
void write(OutputArchiveHandle &handle) const override
Write the object to one or more catalogs.
std::shared_ptr< afw::math::BoundedField > operator*(double const scale) const override
Return a scaled BoundedField.
std::string getPythonModule() const override
Return the fully-qualified Python module that should be imported to guarantee that its factory is reg...
double evaluate(geom::Point2D const &position) const override
Evaluate the field at the given point.
bool operator==(BoundedField const &rhs) const override
BoundedFields (of the same sublcass) are equal if their bounding boxes and parameters are equal.
std::string getPersistenceName() const override
Return the unique name used to persist this object and look up its factory.
Reports arguments outside the domain of an operation.
Definition: Runtime.h:57
Reports attempts to access elements using an invalid key.
Definition: Runtime.h:151
T end(T... args)
T front(T... args)
def scale(algorithm, min, max=None, frame=None)
Definition: ds9.py:108
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
A base class for image defects.
T push_back(T... args)
T reserve(T... args)
T size(T... args)
afw::table::Key< double > weight
afw::table::PointKey< int > bboxMax
afw::table::Key< afw::table::Flag > throwOnMissing
afw::table::Key< int > coaddWcs
afw::table::Key< int > field
afw::table::Key< double > default_
afw::table::Key< int > validPolygon
afw::table::Schema schema
afw::table::PointKey< int > bboxMin
afw::table::Key< int > wcs
Struct used to hold one Exposure's data in a CoaddBoundedField.
std::shared_ptr< afw::geom::SkyWcs const > wcs
std::shared_ptr< afw::math::BoundedField > field
std::shared_ptr< afw::geom::polygon::Polygon const > validPolygon