LSST Applications g063fba187b+fee0456c91,g0f08755f38+ea96e5a5a3,g1653933729+a8ce1bb630,g168dd56ebc+a8ce1bb630,g1a2382251a+90257ff92a,g20f6ffc8e0+ea96e5a5a3,g217e2c1bcf+937a289c59,g28da252d5a+daa7da44eb,g2bbee38e9b+253935c60e,g2bc492864f+253935c60e,g3156d2b45e+6e55a43351,g32e5bea42b+31359a2a7a,g347aa1857d+253935c60e,g35bb328faa+a8ce1bb630,g3a166c0a6a+253935c60e,g3b1af351f3+a8ce1bb630,g3e281a1b8c+c5dd892a6c,g414038480c+416496e02f,g41af890bb2+afe91b1188,g599934f4f4+0db33f7991,g7af13505b9+e36de7bce6,g80478fca09+da231ba887,g82479be7b0+a4516e59e3,g858d7b2824+ea96e5a5a3,g89c8672015+f4add4ffd5,g9125e01d80+a8ce1bb630,ga5288a1d22+bc6ab8dfbd,gb58c049af0+d64f4d3760,gc28159a63d+253935c60e,gcab2d0539d+3f2b72788c,gcf0d15dbbd+4ea9c45075,gda6a2b7d83+4ea9c45075,gdaeeff99f8+1711a396fd,ge79ae78c31+253935c60e,gef2f8181fd+3031e3cf99,gf0baf85859+c1f95f4921,gfa517265be+ea96e5a5a3,gfa999e8aa5+17cd334064,w.2024.50
LSST Data Management Base Package
Loading...
Searching...
No Matches
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
31namespace lsst {
32namespace afw {
33namespace table {
34namespace io {
35
39 );
40
41} // namespace io
42} // namespace table
43
44namespace math {
45
46PixelAreaBoundedField::PixelAreaBoundedField(
49 lsst::geom::AngleUnit const & unit,
50 double scaling
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
68ndarray::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
121 os << "PixelAreaBoundedField(" << (*_skyWcs) << ", scaling=" << _scaling << ")";
122 return os.str();
123}
124
125
126namespace {
127
128struct 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
139private:
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
154class PixelAreaBoundedFieldFactory : public table::io::PersistableFactory {
155public:
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
173std::string getPixelAreaBoundedFieldPersistenceName() { return "PixelAreaBoundedField"; }
174
175PixelAreaBoundedFieldFactory registration(getPixelAreaBoundedFieldPersistenceName());
176
177} // namespace
178
180 return getPixelAreaBoundedFieldPersistenceName();
181}
182
183std::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
AmpInfoBoxKey bbox
Definition Amplifier.cc:117
#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
A BoundedField that gives the amount a pixel is distorted at each point.
std::ostream * os
Definition Schema.cc:557
table::Key< table::Array< std::uint8_t > > wcs
Definition SkyWcs.cc:66
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
table::Schema schema
Definition python.h:134
An abstract base class for 2-d functions defined on an integer bounding boxes.
lsst::geom::Box2I getBBox() const
Return the bounding box that defines the region where the field is valid.
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.
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.
A class used to convert scalar POD types such as double to Angle.
Definition Angle.h:71
An integer coordinate rectangle.
Definition Box.h:55
Reports invalid arguments.
Definition Runtime.h:66
T emplace_back(T... args)
BoxKey< lsst::geom::Box2I > Box2IKey
Definition aggregates.h:283
AngleUnit constexpr radians
constant with units of radians
Definition Angle.h:109
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