LSST Applications g04e9c324dd+8c5ae1fdc5,g134cb467dc+1b3060144d,g18429d2f64+f642bf4753,g199a45376c+0ba108daf9,g1fd858c14a+2dcf163641,g262e1987ae+7b8c96d2ca,g29ae962dfc+3bd6ecb08a,g2cef7863aa+aef1011c0b,g35bb328faa+8c5ae1fdc5,g3fd5ace14f+53e1a9e7c5,g4595892280+fef73a337f,g47891489e3+2efcf17695,g4d44eb3520+642b70b07e,g53246c7159+8c5ae1fdc5,g67b6fd64d1+2efcf17695,g67fd3c3899+b70e05ef52,g74acd417e5+317eb4c7d4,g786e29fd12+668abc6043,g87389fa792+8856018cbb,g89139ef638+2efcf17695,g8d7436a09f+3be3c13596,g8ea07a8fe4+9f5ccc88ac,g90f42f885a+a4e7b16d9b,g97be763408+ad77d7208f,g9dd6db0277+b70e05ef52,ga681d05dcb+a3f46e7fff,gabf8522325+735880ea63,gac2eed3f23+2efcf17695,gb89ab40317+2efcf17695,gbf99507273+8c5ae1fdc5,gd8ff7fe66e+b70e05ef52,gdab6d2f7ff+317eb4c7d4,gdc713202bf+b70e05ef52,gdfd2d52018+b10e285e0f,ge365c994fd+310e8507c4,ge410e46f29+2efcf17695,geaed405ab2+562b3308c0,gffca2db377+8c5ae1fdc5,w.2025.35
LSST Data Management Base Package
Loading...
Searching...
No Matches
CentroidUtilities.cc
Go to the documentation of this file.
1/*
2 * LSST Data Management System
3 * Copyright 2008-2014 LSST Corporation.
4 *
5 * This product includes software developed by the
6 * LSST Project (http://www.lsst.org/).
7 *
8 * This program is free software: you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation, either version 3 of the License, or
11 * (at your option) any later version.
12 *
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the LSST License Statement and
19 * the GNU General Public License along with this program. If not,
20 * see <http://www.lsstcorp.org/LegalNotices/>.
21 */
22
23#include <cmath>
24
25#include "lsst/geom/Angle.h"
26#include "lsst/geom/Point.h"
29
30namespace lsst {
31namespace meas {
32namespace base {
33
35 : x(std::numeric_limits<CentroidElement>::quiet_NaN()),
36 y(std::numeric_limits<CentroidElement>::quiet_NaN()),
37 xErr(std::numeric_limits<ErrElement>::quiet_NaN()),
38 yErr(std::numeric_limits<ErrElement>::quiet_NaN()),
39 x_y_Cov(std::numeric_limits<ErrElement>::quiet_NaN()) {}
40
41Centroid const CentroidResult::getCentroid() const { return Centroid(x, y); }
42
44 x = centroid.getX();
45 y = centroid.getY();
46}
47
50 m << xErr * xErr, x_y_Cov, x_y_Cov, yErr * yErr;
51 return m;
52}
53
55 xErr = std::sqrt(matrix(0, 0));
56 yErr = std::sqrt(matrix(1, 1));
57 x_y_Cov = matrix(0, 1);
58}
59
61 xErr = _xErr;
62 yErr = _yErr;
63 x_y_Cov = 0.0;
64}
65
67 os << "x=" << result.x << ", y=" << result.y << ", xErr=" << result.xErr << ", yErr=" << result.yErr;
68 return os;
69}
70
72 std::string const &doc, UncertaintyEnum uncertainty) {
74 r._centroid = afw::table::PointKey<CentroidElement>::addFields(schema, name, doc, "pixel");
75 if (uncertainty != NO_UNCERTAINTY) {
78 sigma[0] = schema.addField<ErrElement>(schema.join(name, "xErr"), "1-sigma uncertainty on x position",
79 "pixel");
80 sigma[1] = schema.addField<ErrElement>(schema.join(name, "yErr"), "1-sigma uncertainty on y position",
81 "pixel");
82 if (uncertainty == FULL_COVARIANCE) {
83 cov.push_back(schema.addField<ErrElement>(schema.join(name, "x_y_Cov"),
84 "uncertainty covariance in x and y", "pixel^2"));
85 }
86 r._centroidErr = afw::table::CovarianceMatrixKey<ErrElement, 2>(sigma, cov);
87 }
88 return r;
89}
90
91namespace {
92
93std::vector<std::string> getNameVector() {
95 v.push_back("x");
96 v.push_back("y");
97 return v;
98}
99
100} // namespace
101
103 static std::vector<std::string> names = getNameVector(); // C++11 TODO: just use initializer list
104 try {
105 _centroidErr = afw::table::CovarianceMatrixKey<ErrElement, 2>(s, names);
107 }
108}
109
112 r.setCentroid(record.get(_centroid));
113 if (_centroidErr.isValid()) {
114 r.setCentroidErr(record.get(_centroidErr));
115 }
116 return r;
117}
118
120 record.set(_centroid, value.getCentroid());
121 if (_centroidErr.isValid()) {
122 record.set(_centroidErr, value.getCentroidErr());
123 }
124}
125
127 : BaseTransform{name} {
128 // Map the flag through to the output
129 mapper.addMapping(mapper.getInputSchema().find<afw::table::Flag>(name + "_flag").key);
130
131 // Add keys for the coordinates
132 auto &s = mapper.editOutputSchema();
133 _coordKey = afw::table::CoordKey::addFields(s, name, "ICRS coordinates");
134
135 // If the centroid has an error key we also include one on the celestial
136 // coordinates; otherwise, it isn't necessary. Note that if we provide for
137 // errors in celestial coordinates, we always need the full covariance.
138 if (CentroidResultKey(mapper.getInputSchema()[name]).getCentroidErr().isValid()) {
141 sigma[0] = s.addField<ErrElement>(s.join(name, "raErr"), "1-sigma uncertainty on RA", "rad");
142 sigma[1] = s.addField<ErrElement>(s.join(name, "decErr"), "1-sigma uncertainty on dec", "rad");
143 cov[0] = s.addField<ErrElement>(s.join(name, "ra_dec_Cov"), "Uncertainty covariance in RA and dec",
144 "rad^2");
145 _coordErrKey = afw::table::CovarianceMatrixKey<ErrElement, 2>(sigma, cov);
146 }
147}
148
150 afw::table::BaseCatalog &outputCatalog, afw::geom::SkyWcs const &wcs,
151 afw::image::PhotoCalib const &photoCalib) const {
152 checkCatalogSize(inputCatalog, outputCatalog);
153 CentroidResultKey centroidResultKey(inputCatalog.getSchema()[_name]);
154
155 afw::table::SourceCatalog::const_iterator inSrc = inputCatalog.begin();
156 afw::table::BaseCatalog::iterator outSrc = outputCatalog.begin();
157
158 for (; inSrc != inputCatalog.end() && outSrc != outputCatalog.end(); ++inSrc, ++outSrc) {
159 CentroidResult centroidResult = centroidResultKey.get(*inSrc);
160
161 _coordKey.set(*outSrc, wcs.pixelToSky(centroidResult.getCentroid()));
162
163 if (centroidResultKey.getCentroidErr().isValid()) {
164 CentroidCov centroidCov = centroidResult.getCentroidErr();
165 if (!(std::isnan(centroidCov(0, 0)) || std::isnan(centroidCov(1, 1)))) {
166 auto transform = wcs.linearizePixelToSky(centroidResult.getCentroid(), geom::radians)
167 .getLinear()
168 .getMatrix();
169 _coordErrKey.set(*outSrc, (transform * centroidResult.getCentroidErr().cast<double>() *
170 transform.transpose())
171 .cast<ErrElement>());
172 }
173 }
174 }
175}
176
177// Add a key "flag_resetToPeak" if something is wrong with the centroid
178// Save a key to this algorithm's Centroid, as well as the general failure flag
179CentroidChecker::CentroidChecker(afw::table::Schema &schema, std::string const &name, bool doFootprintCheck,
180 double maxDistFromPeak)
181 : _doFootprintCheck(doFootprintCheck), _maxDistFromPeak(maxDistFromPeak) {
182 _resetKey = schema.addField<afw::table::Flag>(schema.join(name, "flag_resetToPeak"),
183 "set if CentroidChecker reset the centroid");
184 _failureKey = schema.find<afw::table::Flag>(schema.join(name, "flag")).key;
185 _xKey = schema.find<CentroidElement>(schema.join(name, "x")).key;
186 _yKey = schema.find<CentroidElement>(schema.join(name, "y")).key;
187
188 // We only check errors on the centroid if they exist: not all measurement
189 // algorithms provide them.
190 try {
191 _xErrKey = schema.find<ErrElement>(schema.join(name, "xErr")).key;
192 } catch (pex::exceptions::NotFoundError &err) {
193 }
194 try {
195 _yErrKey = schema.find<ErrElement>(schema.join(name, "yErr")).key;
196 } catch (pex::exceptions::NotFoundError &err) {
197 }
198 if (_xErrKey.isValid() || _yErrKey.isValid()) {
199 _badErrorKey = schema.addField<afw::table::Flag>(schema.join(name, "flag_badError"),
200 "Error on x and/or y position is NaN");
201 }
202}
203
205 CentroidElement x = record.get(_xKey);
206 CentroidElement y = record.get(_yKey);
207
208 // Check any errors specified on the centroid position are valid.
209 if ((_xErrKey.isValid() && std::isnan(record.get(_xErrKey))) ||
210 (_yErrKey.isValid() && std::isnan(record.get(_yErrKey)))) {
211 record.set(_badErrorKey, true);
212 record.set(_failureKey, true);
213 }
214
215 // Only proceed with checking if appropriately configured.
216 if (!_doFootprintCheck && _maxDistFromPeak < 0.0) {
217 return false;
218 }
219
220 // Check that the centroid has a footprint that we can validate; otherwise, give up.
222 if (!footprint) {
223 throw LSST_EXCEPT(pex::exceptions::RuntimeError, "No Footprint attached to record");
224 }
225 if (footprint->getPeaks().empty()) {
226 throw LSST_EXCEPT(pex::exceptions::RuntimeError, "Footprint has no peaks; cannot verify centroid.");
227 }
228
229 // Set the centroid to the first footprint if the centroid is either more than
230 // _maxDistFromPeak pixels from the centroid, or if it is outside the footprint.
231 CentroidElement footX = footprint->getPeaks().front().getFx();
232 CentroidElement footY = footprint->getPeaks().front().getFy();
233 double distsq = (x - footX) * (x - footX) + (y - footY) * (y - footY);
234 if ((_doFootprintCheck && !footprint->contains(geom::Point2I(geom::Point2D(x, y)))) ||
235 ((_maxDistFromPeak > 0) && (distsq > _maxDistFromPeak * _maxDistFromPeak))) {
236 record.set(_xKey, footX);
237 record.set(_yKey, footY);
238 record.set(_failureKey, true);
239 record.set(_resetKey, true);
240 return true;
241 }
242 return false;
243}
244} // namespace base
245} // namespace meas
246} // namespace lsst
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
A 2-dimensional celestial WCS that transform pixels to ICRS RA/Dec, using the LSST standard for pixel...
Definition SkyWcs.h:118
The photometric calibration of an exposure.
Definition PhotoCalib.h:114
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
void set(Key< T > const &key, U const &value)
Set value of a field for the given key.
Definition BaseRecord.h:164
CatalogIterator< typename Internal::iterator > iterator
Definition Catalog.h:110
iterator begin()
Iterator access.
Definition Catalog.h:401
static CoordKey addFields(afw::table::Schema &schema, std::string const &name, std::string const &doc)
Add a pair of _ra, _dec fields to a Schema, and return a CoordKey that points to them.
bool isValid() const noexcept
Return True if all the constituent error Keys are valid.
static PointKey addFields(Schema &schema, std::string const &name, std::string const &doc, std::string const &unit)
Add a pair of _x, _y fields to a Schema, and return a PointKey that points to them.
Definition aggregates.cc:36
Defines the fields and offsets for a table.
Definition Schema.h:51
SchemaItem< T > find(std::string const &name) const
Find a SchemaItem in the Schema by name.
Definition Schema.cc:467
A mapping between the keys of two Schemas, used to copy data between them.
Key< T > addMapping(Key< T > const &inputKey, bool doReplace=false)
Add a new field to the output Schema that is a copy of a field in the input Schema.
Schema & editOutputSchema()
Return a reference to the output schema that allows it to be modified in place.
Schema const getInputSchema() const
Return the input schema (copy-on-write).
typename Base::const_iterator const_iterator
Record class that contains measurements made on a single exposure.
Definition Source.h:78
std::shared_ptr< Footprint > getFootprint() const
Definition Source.h:98
A proxy type for name lookups in a Schema.
Definition Schema.h:367
void checkCatalogSize(afw::table::BaseCatalog const &cat1, afw::table::BaseCatalog const &cat2) const
Ensure that catalogs have the same size.
Definition Transform.h:102
BaseTransform(std::string const &name)
Definition Transform.h:88
CentroidChecker(afw::table::Schema &schema, std::string const &name, bool inside=true, double maxDistFromPeak=-1.0)
Check source record produced by a centroid algorithm called "name".
bool operator()(afw::table::SourceRecord &record) const
Set the centroid to the first footprint if the centroid is either more than _dist pixels from the foo...
A FunctorKey for CentroidResult.
virtual void set(afw::table::BaseRecord &record, CentroidResult const &value) const
Set a CentroidResult in the given record.
afw::table::CovarianceMatrixKey< ErrElement, 2 > getCentroidErr() const
Return a FunctorKey to just the uncertainty matrix.
CentroidResultKey()
Default constructor; instance will not be usuable unless subsequently assigned to.
virtual CentroidResult get(afw::table::BaseRecord const &record) const
Get a CentroidResult from the given record.
static CentroidResultKey addFields(afw::table::Schema &schema, std::string const &name, std::string const &doc, UncertaintyEnum uncertainty)
Add the appropriate fields to a Schema, and return a CentroidResultKey that manages them.
CentroidTransform(std::string const &name, afw::table::SchemaMapper &mapper)
virtual void operator()(afw::table::SourceCatalog const &inputCatalog, afw::table::BaseCatalog &outputCatalog, afw::geom::SkyWcs const &wcs, afw::image::PhotoCalib const &photoCalib) const
Reports attempts to access elements using an invalid key.
Definition Runtime.h:151
Reports errors that are due to events beyond the control of the program.
Definition Runtime.h:104
T isnan(T... args)
CatalogT< BaseRecord > BaseCatalog
Definition fwd.h:72
AngleUnit constexpr radians
constant with units of radians
Definition Angle.h:109
Point< double, 2 > Point2D
Definition Point.h:324
Point< int, 2 > Point2I
Definition Point.h:321
UncertaintyEnum
An enum used to specify how much uncertainty information measurement algorithms provide.
Definition constants.h:43
@ FULL_COVARIANCE
The full covariance matrix is provided.
Definition constants.h:46
@ NO_UNCERTAINTY
Algorithm provides no uncertainy information at all.
Definition constants.h:44
std::ostream & operator<<(std::ostream &os, CentroidResult const &result)
Eigen::Matrix< ErrElement, 2, 2, Eigen::DontAlign > CentroidCov
Definition constants.h:59
double CentroidElement
Definition constants.h:56
geom::Point< CentroidElement, 2 > Centroid
Definition constants.h:58
STL namespace.
T push_back(T... args)
T sqrt(T... args)
A reusable struct for centroid measurements.
CentroidElement y
y (row) coordinate of the measured position
Centroid const getCentroid() const
Return a Point object containing the measured x and y.
CentroidElement x
x (column) coordinate of the measured position
void setCentroidErr(CentroidCov const &matrix)
Set the struct uncertainty fields from the given matrix, with rows and columns ordered (x,...
CentroidCov const getCentroidErr() const
Return the 2x2 symmetric covariance matrix, with rows and columns ordered (x, y)
void setCentroid(Centroid const &centroid)
Set the struct fields from the given Point object.
ErrElement yErr
standard deviation of y
CentroidResult()
Constructor; initializes everything to NaN.
ErrElement x_y_Cov
x,y term in the uncertainty convariance matrix
ErrElement xErr
standard deviation of x