LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
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 
30 
31 namespace lsst { namespace meas { namespace algorithms {
32 
33 
35  afw::geom::Box2I const & bbox,
37  ElementVector const & elements
38 ) :
39  afw::math::BoundedField(bbox),
40  _throwOnMissing(true),
41  _default(0.0), // unused
42  _coaddWcs(coaddWcs),
43  _elements(elements)
44 {}
45 
47  afw::geom::Box2I const & bbox,
49  ElementVector const & elements,
50  double default_
51 ) :
52  afw::math::BoundedField(bbox),
53  _throwOnMissing(false),
54  _default(default_),
55  _coaddWcs(coaddWcs),
56  _elements(elements)
57 {}
58 
59 double CoaddBoundedField::evaluate(afw::geom::Point2D const & position) const {
60  PTR(afw::coord::Coord) coord = _coaddWcs->pixelToSky(position);
61  double sum = 0.0;
62  double wSum = 0.0;
63  for (ElementVector::const_iterator i = _elements.begin(); i != _elements.end(); ++i) {
64  afw::geom::Point2D transformedPosition = i->wcs->skyToPixel(*coord);
65  bool inValidArea = i->validPolygon ? i->validPolygon->contains(transformedPosition) : true;
66  if (afw::geom::Box2D(i->field->getBBox()).contains(transformedPosition) && inValidArea) {
67  sum += i->weight * i->field->evaluate(transformedPosition);
68  wSum += i->weight;
69  }
70  }
71  if (wSum == 0.0) {
72  if (_throwOnMissing) {
73  throw LSST_EXCEPT(
74  pex::exceptions::DomainError,
75  (boost::format("No constituent fields to evaluate at point %f, %f")
76  % position.getX() % position.getY()).str()
77  );
78  } else {
79  return _default;
80  }
81  }
82  return sum / wSum;
83 }
84 
85 // ---------- Persistence -----------------------------------------------------------------------------------
86 
87 // For persistence of CoaddBoundedField, we have two catalogs: the first has just one record, and contains
88 // the archive ID of the coadd WCS and the parameters that control missing data. The other catalog
89 // has one record for each element, with fields corresponding to the data members of the Element struct.
90 
91 namespace {
92 
93 namespace tbl = afw::table;
94 
95 // Singleton class that manages the first persistence catalog's schema and keys
96 class CoaddBoundedFieldPersistenceKeys1 : private boost::noncopyable {
97 public:
98  tbl::Schema schema;
99  tbl::PointKey<int> bboxMin;
100  tbl::PointKey<int> bboxMax;
101  tbl::Key<int> coaddWcs;
102  tbl::Key<tbl::Flag> throwOnMissing;
103  tbl::Key<double > default_;
104 
105  static CoaddBoundedFieldPersistenceKeys1 const & get() {
106  static CoaddBoundedFieldPersistenceKeys1 const instance;
107  return instance;
108  }
109 
110 private:
111  CoaddBoundedFieldPersistenceKeys1() :
112  schema(),
113  bboxMin(tbl::PointKey<int>::addFields(
114  schema, "bbox_min", "lower-left corner of bounding box", "pixels")),
115  bboxMax(tbl::PointKey<int>::addFields(
116  schema, "bbox_max", "upper-right corner of bounding box", "pixels")),
117  coaddWcs(schema.addField<int>("coaddWcs", "archive ID of the coadd's WCS")),
118  throwOnMissing(schema.addField<tbl::Flag>("throwOnMissing",
119  "whether to throw an exception on missing data")),
120  default_(schema.addField<double>("default", "default value to use when throwOnMissing is False"))
121  {
122  schema.getCitizen().markPersistent();
123  }
124 };
125 
126 // Singleton class that manages the second persistence catalog's schema and keys
127 class CoaddBoundedFieldPersistenceKeys2 : private boost::noncopyable {
128 public:
129  tbl::Schema schema;
130  tbl::Key<int> field;
131  tbl::Key<int> wcs;
132  tbl::Key<int> validPolygon;
133  tbl::Key<double> weight;
134 
135  static CoaddBoundedFieldPersistenceKeys2 const & get() {
136  static CoaddBoundedFieldPersistenceKeys2 const instance;
137  return instance;
138  }
139 
140 private:
141  CoaddBoundedFieldPersistenceKeys2() :
142  schema(),
143  field(schema.addField<int>("field", "archive ID of the BoundedField to be coadded")),
144  wcs(schema.addField<int>("wcs", "archive ID of the Wcs associated with this element")),
145  validPolygon(schema.addField<int>("validPolygon", "archive ID of the Polygon associated with this element")),
146  weight(schema.addField<double>("weight", "weight value for this element"))
147  {
148  schema.getCitizen().markPersistent();
149  }
150 };
151 
152 } // anonymous
153 
154 class CoaddBoundedField::Factory : public tbl::io::PersistableFactory {
155 public:
156 
157  virtual PTR(tbl::io::Persistable)
158  read(InputArchive const & archive, CatalogVector const & catalogs) const {
159  CoaddBoundedFieldPersistenceKeys1 const & keys1 = CoaddBoundedFieldPersistenceKeys1::get();
160  CoaddBoundedFieldPersistenceKeys2 const & keys2 = CoaddBoundedFieldPersistenceKeys2::get();
161  LSST_ARCHIVE_ASSERT(catalogs.size() == 2u);
162  LSST_ARCHIVE_ASSERT(catalogs.front().getSchema() == keys1.schema);
163  LSST_ARCHIVE_ASSERT(catalogs.back().getSchema() == keys2.schema);
164  tbl::BaseRecord const & record1 = catalogs.front().front();
165  ElementVector elements;
166  elements.reserve(catalogs.back().size());
167  for (tbl::BaseCatalog::const_iterator i = catalogs.back().begin(); i != catalogs.back().end(); ++i) {
168  elements.push_back(
169  Element(
170  archive.get<afw::math::BoundedField>(i->get(keys2.field)),
171  archive.get<afw::image::Wcs>(i->get(keys2.wcs)),
172  archive.get<afw::geom::polygon::Polygon>(i->get(keys2.validPolygon)),
173  i->get(keys2.weight)
174  )
175  );
176  }
177  return boost::make_shared<CoaddBoundedField>(
178  afw::geom::Box2I(record1.get(keys1.bboxMin), record1.get(keys1.bboxMax)),
179  archive.get<afw::image::Wcs>(record1.get(keys1.coaddWcs)),
180  elements,
181  record1.get(keys1.default_)
182  );
183  }
184 
185  Factory(std::string const & name) : tbl::io::PersistableFactory(name) {}
186 
187 };
188 
189 namespace {
190 
191 std::string getCoaddBoundedFieldPersistenceName() { return "CoaddBoundedField"; }
192 
193 CoaddBoundedField::Factory registration(getCoaddBoundedFieldPersistenceName());
194 
195 } // anonymous
196 
197 std::string CoaddBoundedField::getPersistenceName() const { return getCoaddBoundedFieldPersistenceName(); }
198 
199 std::string CoaddBoundedField::getPythonModule() const { return "lsst.meas.algorithms"; }
200 
202  CoaddBoundedFieldPersistenceKeys1 const & keys1 = CoaddBoundedFieldPersistenceKeys1::get();
203  CoaddBoundedFieldPersistenceKeys2 const & keys2 = CoaddBoundedFieldPersistenceKeys2::get();
204  tbl::BaseCatalog cat1 = handle.makeCatalog(keys1.schema);
205  PTR(tbl::BaseRecord) record1 = cat1.addNew();
206  record1->set(keys1.bboxMin, getBBox().getMin());
207  record1->set(keys1.bboxMax, getBBox().getMax());
208  record1->set(keys1.coaddWcs, handle.put(_coaddWcs));
209  record1->set(keys1.default_, _default);
210  handle.saveCatalog(cat1);
211  tbl::BaseCatalog cat2 = handle.makeCatalog(keys2.schema);
212  for (ElementVector::const_iterator i = _elements.begin(); i != _elements.end(); ++i) {
213  PTR(tbl::BaseRecord) record2 = cat2.addNew();
214  record2->set(keys2.field, handle.put(i->field));
215  record2->set(keys2.wcs, handle.put(i->wcs));
216  record2->set(keys2.validPolygon, handle.put(i->validPolygon));
217  record2->set(keys2.weight, i->weight);
218  }
219  handle.saveCatalog(cat2);
220 }
221 
222 PTR(afw::math::BoundedField) CoaddBoundedField::operator*(double const scale) const {
223  throw LSST_EXCEPT(pex::exceptions::NotFoundError, "Scaling of CoaddBoundedField is not implemented");
224 }
225 
226 }}} // namespace lsst::meas::algorithms
227 
228 
#define PTR(...)
Definition: base.h:41
table::Key< std::string > name
Definition: ApCorrMap.cc:71
BaseCatalog makeCatalog(Schema const &schema)
Return a new, empty catalog with the given schema.
boost::shared_ptr< afw::image::Wcs const > _coaddWcs
tbl::Key< double > weight
An object passed to Persistable::write to allow it to persist itself.
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
Definition: Persistable.h:47
int put(Persistable const *obj, bool permissive=false)
Save a nested Persistable to the same archive.
A custom container class for records, based on std::vector.
Definition: Catalog.h:94
virtual void write(OutputArchiveHandle &handle) const
Write the object to one or more catalogs.
virtual boost::shared_ptr< tbl::io::Persistable > read(InputArchive const &archive, CatalogVector const &catalogs) const
virtual double evaluate(afw::geom::Point2D const &position) const
tbl::Key< int > wcs
virtual std::string getPythonModule() const
Return the fully-qualified Python module that should be imported to guarantee that its factory is reg...
Implementation of the WCS standard for a any projection.
Definition: Wcs.h:107
boost::enable_if< typename ExpressionTraits< Scalar >::IsScalar, Scalar >::type sum(Scalar const &scalar)
Definition: operators.h:1250
boost::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:469
An integer coordinate rectangle.
Definition: Box.h:53
tbl::Key< int > validPolygon
CoaddBoundedField(afw::geom::Box2I const &bbox, boost::shared_ptr< afw::image::Wcs const > coaddWcs, ElementVector const &elements)
tbl::Key< int > field
tbl::Key< double > default_
geom::Box2I getBBox() const
Definition: BoundedField.h:95
virtual std::string getPersistenceName() const
Return the unique name used to persist this object and look up its factory.
tbl::Schema schema
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
tbl::PointKey< int > bboxMax
Base class for all records.
Definition: BaseRecord.h:27
tbl::Key< tbl::Flag > throwOnMissing
An abstract base class for 2-d functions defined on an integer bounding boxes.
Definition: BoundedField.h:49
A floating-point coordinate rectangle geometry.
Definition: Box.h:271
void saveCatalog(BaseCatalog const &catalog)
Save a catalog in the archive.
tbl::Key< int > coaddWcs
tbl::PointKey< int > bboxMin