LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
PixelAreaBoundedField.cc
Go to the documentation of this file.
1 /*
2  * This file is part of afw.
3  *
4  * Developed for the LSST Data Management System.
5  * This product includes software developed by the LSST Project
6  * (https://www.lsst.org).
7  * See the COPYRIGHT file at the top-level directory of this distribution
8  * for details of code ownership.
9  *
10  * This program is free software: you can redistribute it and/or modify
11  * it under the terms of the GNU General Public License as published by
12  * the Free Software Foundation, either version 3 of the License, or
13  * (at your option) any later version.
14  *
15  * This program is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  * GNU General Public License for more details.
19  *
20  * You should have received a copy of the GNU General Public License
21  * along with this program. If not, see <https://www.gnu.org/licenses/>.
22  */
23 
28 
30 
31 namespace lsst {
32 namespace afw {
33 namespace table {
34 namespace io {
35 
39  );
40 
41 } // namespace io
42 } // namespace table
43 
44 namespace math {
45 
46 PixelAreaBoundedField::PixelAreaBoundedField(
47  lsst::geom::Box2I const &bbox,
49  lsst::geom::AngleUnit const & unit,
50  double scaling
51 ) : BoundedField(bbox),
52  _skyWcs(skyWcs),
53  _scaling(scaling)
54 {
55  if (_skyWcs == nullptr) {
56  throw LSST_EXCEPT(
58  "SkyWcs passed to PixelAreaBoundedField is null."
59  );
60  }
61  _scaling /= std::pow((1.0*unit).asRadians(), 2);
62 }
63 
65  return std::pow(_skyWcs->getPixelScale(position).asRadians(), 2) * _scaling;
66 }
67 
68 ndarray::Array<double, 1, 1> PixelAreaBoundedField::evaluate(
69  ndarray::Array<double const, 1> const & x,
70  ndarray::Array<double const, 1> const & y
71 ) const {
72  if (x.getShape() != y.getShape()) {
73  throw LSST_EXCEPT(
75  (boost::format("Inconsistent shapes in evaluate; %s != %s.") % x.getShape() % y.getShape()).str()
76  );
77  }
78  // Compute _skyWcs->pixelToSky at all of the given points in a single
79  // vectorized call, along with points one pixel away in x and y.
80  double constexpr side = 1.0;
81  std::size_t const n = x.size();
83  pixPoints.reserve(n*3);
84  for (std::size_t i = 0; i < n; ++i) {
85  pixPoints.emplace_back(x[i], y[i]);
86  pixPoints.emplace_back(x[i] + side, y[i]);
87  pixPoints.emplace_back(x[i], y[i] + side);
88  }
89  auto skyPoints = _skyWcs->pixelToSky(pixPoints);
90  // Work in 3-space to avoid RA wrapping and pole issues.
91  ndarray::Array<double, 1, 1> z = ndarray::allocate(x.getShape());
92  for (std::size_t i = 0; i < n; ++i) {
93  std::size_t j = i*3;
94  auto skyLL = skyPoints[j].getVector();
95  auto skyDx = skyPoints[j + 1].getVector() - skyLL;
96  auto skyDy = skyPoints[j + 2].getVector() - skyLL;
97  double skyAreaSq = skyDx.cross(skyDy).getSquaredNorm();
98  z[i] = _scaling * std::sqrt(skyAreaSq) / (side*side);
99  }
100  return z;
101 }
102 
104  return _skyWcs->isPersistable();
105 }
106 
108  return std::make_shared<PixelAreaBoundedField>(getBBox(), _skyWcs, lsst::geom::radians, _scaling*scale);
109 }
110 
112  auto rhsCasted = dynamic_cast<PixelAreaBoundedField const *>(&rhs);
113  if (!rhsCasted) return false;
114 
115  return getBBox() == rhsCasted->getBBox() && *_skyWcs == *rhsCasted->_skyWcs &&
116  _scaling == rhsCasted->_scaling;
117 }
118 
119 std::string PixelAreaBoundedField::toString() const {
121  os << "PixelAreaBoundedField(" << (*_skyWcs) << ", scaling=" << _scaling << ")";
122  return os.str();
123 }
124 
125 
126 namespace {
127 
128 struct PersistenceHelper {
129  table::Schema schema;
131  table::Key<int> wcs;
132  table::Key<double> scaling;
133 
134  static PersistenceHelper const & get() {
135  static PersistenceHelper const instance;
136  return instance;
137  }
138 
139 private:
140  PersistenceHelper() :
141  schema(),
142  bbox(table::Box2IKey::addFields(schema, "bbox", "Bounding box for field.", "pixel")),
143  wcs(schema.addField<int>("wcs", "Archive ID for SkyWcs instance.")),
144  scaling(schema.addField<double>("scaling",
145  "Scaling factor (including any transformation from rad^2."))
146  {}
147  PersistenceHelper(PersistenceHelper const &) = delete;
148  PersistenceHelper(PersistenceHelper &&) = delete;
149  PersistenceHelper & operator=(PersistenceHelper const &) = delete;
150  PersistenceHelper & operator=(PersistenceHelper &&) = delete;
151  ~PersistenceHelper() noexcept = default;
152 };
153 
154 class PixelAreaBoundedFieldFactory : public table::io::PersistableFactory {
155 public:
156  explicit PixelAreaBoundedFieldFactory(std::string const& name)
157  : afw::table::io::PersistableFactory(name) {}
158 
159  std::shared_ptr<table::io::Persistable> read(InputArchive const& archive,
160  CatalogVector const& catalogs) const override {
161  LSST_ARCHIVE_ASSERT(catalogs.size() == 1u);
162  LSST_ARCHIVE_ASSERT(catalogs.front().size() == 1u);
163  table::BaseRecord const& record = catalogs.front().front();
164  auto const & keys = PersistenceHelper::get();
165  LSST_ARCHIVE_ASSERT(record.getSchema() == keys.schema);
166  lsst::geom::Box2I bbox(record.get(keys.bbox));
167  auto wcs = archive.get<afw::geom::SkyWcs>(record.get(keys.wcs));
168  double scaling = record.get(keys.scaling);
169  return std::make_shared<PixelAreaBoundedField>(bbox, wcs, lsst::geom::radians, scaling);
170  }
171 };
172 
173 std::string getPixelAreaBoundedFieldPersistenceName() { return "PixelAreaBoundedField"; }
174 
175 PixelAreaBoundedFieldFactory registration(getPixelAreaBoundedFieldPersistenceName());
176 
177 } // namespace
178 
180  return getPixelAreaBoundedFieldPersistenceName();
181 }
182 
183 std::string PixelAreaBoundedField::getPythonModule() const { return "lsst.afw.math"; }
184 
186  auto const & keys = PersistenceHelper::get();
187  table::BaseCatalog catalog = handle.makeCatalog(keys.schema);
188  std::shared_ptr<table::BaseRecord> record = catalog.addNew();
189  record->set(keys.bbox, getBBox());
190  record->set(keys.wcs, handle.put(_skyWcs));
191  record->set(keys.scaling, _scaling);
192  handle.saveCatalog(catalog);
193 }
194 
195 } // namespace math
196 } // namespace afw
197 } // namespace lsst
table::Key< std::string > name
Definition: Amplifier.cc:116
double x
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
double z
Definition: Match.cc:44
table::Key< double > scaling
table::Box2IKey bbox
table::Schema schema
table::Key< int > wcs
A BoundedField that gives the amount a pixel is distorted at each point.
std::ostream * os
Definition: Schema.cc:557
int y
Definition: SpanSet.cc:48
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
Definition: Persistable.h:48
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
A BoundedField that evaluate the pixel area of a SkyWcs in angular units.
double evaluate(lsst::geom::Point2D const &position) const override
Evaluate the field at the given point.
std::shared_ptr< BoundedField > operator*(double const scale) const override
Return a scaled BoundedField.
bool isPersistable() const noexcept override
PixelAreaBoundedField is persistable if and only if the nested SkyWcs is.
void write(OutputArchiveHandle &handle) const override
Write the object to one or more catalogs.
std::string getPersistenceName() const override
Return the unique name used to persist this object and look up its factory.
std::string getPythonModule() const override
Return the fully-qualified Python module that should be imported to guarantee that its factory is reg...
bool operator==(BoundedField const &rhs) const override
BoundedFields (of the same sublcass) are equal if their bounding boxes and parameters are equal.
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
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 class used to convert scalar POD types such as double to Angle.
Definition: Angle.h:70
An integer coordinate rectangle.
Definition: Box.h:55
Reports invalid arguments.
Definition: Runtime.h:66
T emplace_back(T... args)
def scale(algorithm, min, max=None, frame=None)
Definition: ds9.py:108
FilterProperty & operator=(FilterProperty const &)=default
BoxKey< lsst::geom::Box2I > Box2IKey
Definition: aggregates.h:201
constexpr AngleUnit radians
constant with units of radians
Definition: Angle.h:108
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
A base class for image defects.
T pow(T... args)
T reserve(T... args)
T sqrt(T... args)
std::shared_ptr< table::io::Persistable > read(table::io::InputArchive const &archive, table::io::CatalogVector const &catalogs) const override
Definition: warpExposure.cc:0