LSST Applications  21.0.0-147-g0e635eb1+1acddb5be5,22.0.0+052faf71bd,22.0.0+1ea9a8b2b2,22.0.0+6312710a6c,22.0.0+729191ecac,22.0.0+7589c3a021,22.0.0+9f079a9461,22.0.1-1-g7d6de66+b8044ec9de,22.0.1-1-g87000a6+536b1ee016,22.0.1-1-g8e32f31+6312710a6c,22.0.1-10-gd060f87+016f7cdc03,22.0.1-12-g9c3108e+df145f6f68,22.0.1-16-g314fa6d+c825727ab8,22.0.1-19-g93a5c75+d23f2fb6d8,22.0.1-19-gb93eaa13+aab3ef7709,22.0.1-2-g8ef0a89+b8044ec9de,22.0.1-2-g92698f7+9f079a9461,22.0.1-2-ga9b0f51+052faf71bd,22.0.1-2-gac51dbf+052faf71bd,22.0.1-2-gb66926d+6312710a6c,22.0.1-2-gcb770ba+09e3807989,22.0.1-20-g32debb5+b8044ec9de,22.0.1-23-gc2439a9a+fb0756638e,22.0.1-3-g496fd5d+09117f784f,22.0.1-3-g59f966b+1e6ba2c031,22.0.1-3-g849a1b8+f8b568069f,22.0.1-3-gaaec9c0+c5c846a8b1,22.0.1-32-g5ddfab5d3+60ce4897b0,22.0.1-4-g037fbe1+64e601228d,22.0.1-4-g8623105+b8044ec9de,22.0.1-5-g096abc9+d18c45d440,22.0.1-5-g15c806e+57f5c03693,22.0.1-7-gba73697+57f5c03693,master-g6e05de7fdc+c1283a92b8,master-g72cdda8301+729191ecac,w.2021.39
LSST Data Management Base Package
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 
30 namespace lsst {
31 namespace meas {
32 namespace 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 
41 Centroid const CentroidResult::getCentroid() const { return Centroid(x, y); }
42 
44  x = centroid.getX();
45  y = centroid.getY();
46 }
47 
49  CentroidCov m;
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  std::string const &doc, UncertaintyEnum uncertainty) {
69  r._centroid = afw::table::PointKey<CentroidElement>::addFields(schema, name, doc, "pixel");
70  if (uncertainty != NO_UNCERTAINTY) {
73  sigma[0] = schema.addField<ErrElement>(schema.join(name, "xErr"), "1-sigma uncertainty on x position",
74  "pixel");
75  sigma[1] = schema.addField<ErrElement>(schema.join(name, "yErr"), "1-sigma uncertainty on y position",
76  "pixel");
77  if (uncertainty == FULL_COVARIANCE) {
78  cov.push_back(schema.addField<ErrElement>(schema.join(name, "x_y_Cov"),
79  "uncertainty covariance in x and y", "pixel^2"));
80  }
82  }
83  return r;
84 }
85 
86 namespace {
87 
88 std::vector<std::string> getNameVector() {
90  v.push_back("x");
91  v.push_back("y");
92  return v;
93 }
94 
95 } // namespace
96 
98  static std::vector<std::string> names = getNameVector(); // C++11 TODO: just use initializer list
99  try {
100  _centroidErr = afw::table::CovarianceMatrixKey<ErrElement, 2>(s, names);
101  } catch (pex::exceptions::NotFoundError &) {
102  }
103 }
104 
106  CentroidResult r;
107  r.setCentroid(record.get(_centroid));
108  if (_centroidErr.isValid()) {
109  r.setCentroidErr(record.get(_centroidErr));
110  }
111  return r;
112 }
113 
115  record.set(_centroid, value.getCentroid());
116  if (_centroidErr.isValid()) {
117  record.set(_centroidErr, value.getCentroidErr());
118  }
119 }
120 
122  : BaseTransform{name} {
123  // Map the flag through to the output
124  mapper.addMapping(mapper.getInputSchema().find<afw::table::Flag>(name + "_flag").key);
125 
126  // Add keys for the coordinates
127  auto &s = mapper.editOutputSchema();
128  _coordKey = afw::table::CoordKey::addFields(s, name, "ICRS coordinates");
129 
130  // If the centroid has an error key we also include one on the celestial
131  // coordinates; otherwise, it isn't necessary. Note that if we provide for
132  // errors in celestial coordinates, we always need the full covariance.
133  if (CentroidResultKey(mapper.getInputSchema()[name]).getCentroidErr().isValid()) {
136  sigma[0] = s.addField<ErrElement>(s.join(name, "raErr"), "1-sigma uncertainty on RA", "rad");
137  sigma[1] = s.addField<ErrElement>(s.join(name, "decErr"), "1-sigma uncertainty on dec", "rad");
138  cov[0] = s.addField<ErrElement>(s.join(name, "ra_dec_Cov"), "Uncertainty covariance in RA and dec",
139  "rad^2");
141  }
142 }
143 
145  afw::table::BaseCatalog &outputCatalog, afw::geom::SkyWcs const &wcs,
146  afw::image::PhotoCalib const &photoCalib) const {
147  checkCatalogSize(inputCatalog, outputCatalog);
148  CentroidResultKey centroidResultKey(inputCatalog.getSchema()[_name]);
149 
150  afw::table::SourceCatalog::const_iterator inSrc = inputCatalog.begin();
151  afw::table::BaseCatalog::iterator outSrc = outputCatalog.begin();
152 
153  for (; inSrc != inputCatalog.end() && outSrc != outputCatalog.end(); ++inSrc, ++outSrc) {
154  CentroidResult centroidResult = centroidResultKey.get(*inSrc);
155 
156  _coordKey.set(*outSrc, wcs.pixelToSky(centroidResult.getCentroid()));
157 
158  if (centroidResultKey.getCentroidErr().isValid()) {
159  CentroidCov centroidCov = centroidResult.getCentroidErr();
160  if (!(std::isnan(centroidCov(0, 0)) || std::isnan(centroidCov(1, 1)))) {
161  auto transform = wcs.linearizePixelToSky(centroidResult.getCentroid(), geom::radians)
162  .getLinear()
163  .getMatrix();
164  _coordErrKey.set(*outSrc, (transform * centroidResult.getCentroidErr().cast<double>() *
165  transform.transpose())
166  .cast<ErrElement>());
167  }
168  }
169  }
170 }
171 
172 // Add a key "flag_resetToPeak" if something is wrong with the centroid
173 // Save a key to this algorithm's Centroid, as well as the general failure flag
175  double maxDistFromPeak)
176  : _doFootprintCheck(doFootprintCheck), _maxDistFromPeak(maxDistFromPeak) {
177  _resetKey = schema.addField<afw::table::Flag>(schema.join(name, "flag_resetToPeak"),
178  "set if CentroidChecker reset the centroid");
179  _failureKey = schema.find<afw::table::Flag>(schema.join(name, "flag")).key;
180  _xKey = schema.find<CentroidElement>(schema.join(name, "x")).key;
181  _yKey = schema.find<CentroidElement>(schema.join(name, "y")).key;
182 
183  // We only check errors on the centroid if they exist: not all measurement
184  // algorithms provide them.
185  try {
186  _xErrKey = schema.find<ErrElement>(schema.join(name, "xErr")).key;
187  } catch (pex::exceptions::NotFoundError &err) {}
188  try {
189  _yErrKey = schema.find<ErrElement>(schema.join(name, "yErr")).key;
190  } catch (pex::exceptions::NotFoundError &err) {}
191  if (_xErrKey.isValid() || _yErrKey.isValid()) {
192  _badErrorKey = schema.addField<afw::table::Flag>(schema.join(name, "flag_badError"),
193  "Error on x and/or y position is NaN");
194  }
195 }
196 
198  CentroidElement x = record.get(_xKey);
199  CentroidElement y = record.get(_yKey);
200 
201  // Check any errors specified on the centroid position are valid.
202  if ((_xErrKey.isValid() && std::isnan(record.get(_xErrKey))) ||
203  (_yErrKey.isValid() && std::isnan(record.get(_yErrKey)))) {
204  record.set(_badErrorKey, true);
205  record.set(_failureKey, true);
206  }
207 
208  // Only proceed with checking if appropriately configured.
209  if (!_doFootprintCheck && _maxDistFromPeak < 0.0) {
210  return false;
211  }
212 
213  // Check that the centroid has a footprint that we can validate; otherwise, give up.
215  if (!footprint) {
216  throw LSST_EXCEPT(pex::exceptions::RuntimeError, "No Footprint attached to record");
217  }
218  if (footprint->getPeaks().empty()) {
219  throw LSST_EXCEPT(pex::exceptions::RuntimeError, "Footprint has no peaks; cannot verify centroid.");
220  }
221 
222  // Set the centroid to the first footprint if the centroid is either more than
223  // _maxDistFromPeak pixels from the centroid, or if it is outside the footprint.
224  CentroidElement footX = footprint->getPeaks().front().getFx();
225  CentroidElement footY = footprint->getPeaks().front().getFy();
226  double distsq = (x - footX) * (x - footX) + (y - footY) * (y - footY);
227  if ((_doFootprintCheck && !footprint->contains(geom::Point2I(geom::Point2D(x, y)))) ||
228  ((_maxDistFromPeak > 0) && (distsq > _maxDistFromPeak * _maxDistFromPeak))) {
229  record.set(_xKey, footX);
230  record.set(_yKey, footY);
231  record.set(_failureKey, true);
232  record.set(_resetKey, true);
233  return true;
234  }
235  return false;
236 }
237 } // namespace base
238 } // namespace meas
239 } // 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
afw::table::Key< double > sigma
Definition: GaussianPsf.cc:49
SchemaMapper * mapper
Definition: SchemaMapper.cc:71
table::Key< table::Array< std::uint8_t > > wcs
Definition: SkyWcs.cc:66
int y
Definition: SpanSet.cc:48
int m
Definition: SpanSet.cc:48
table::Key< int > transform
table::Schema schema
Definition: python.h:134
A 2-dimensional celestial WCS that transform pixels to ICRS RA/Dec, using the LSST standard for pixel...
Definition: SkyWcs.h:117
The photometric calibration of an exposure.
Definition: PhotoCalib.h:114
Base class for all records.
Definition: BaseRecord.h:31
void set(Key< T > const &key, U const &value)
Set value of a field for the given key.
Definition: BaseRecord.h:164
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
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.
Definition: aggregates.cc:83
void set(BaseRecord &record, lsst::geom::SpherePoint const &value) const override
Set an lsst::geom::SpherePoint in the given record.
Definition: aggregates.cc:93
void set(BaseRecord &record, Eigen::Matrix< T, N, N > const &value) const override
Set a covariance matrix in the given record (uses only the lower triangle of the given matrix)
Definition: aggregates.cc:278
bool isValid() const noexcept
Return True if all the constituent error Keys are valid.
Definition: aggregates.cc:294
bool isValid() const noexcept
Return true if the key was initialized to valid offset.
Definition: Key.h:97
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
A mapping between the keys of two Schemas, used to copy data between them.
Definition: SchemaMapper.h:21
typename Base::const_iterator const_iterator
Definition: SortedCatalog.h:50
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
A coordinate class intended to represent absolute positions.
Definition: Point.h:169
Abstract base class for all C++ measurement transformations.
Definition: Transform.h:86
void checkCatalogSize(afw::table::BaseCatalog const &cat1, afw::table::BaseCatalog const &cat2) const
Ensure that catalogs have the same size.
Definition: Transform.h:102
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.
afw::table::CovarianceMatrixKey< ErrElement, 2 > getCentroidErr() const
Return a FunctorKey to just the uncertainty matrix.
virtual void set(afw::table::BaseRecord &record, CentroidResult const &value) const
Set a CentroidResult in the given record.
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)
constexpr AngleUnit radians
constant with units of radians
Definition: Angle.h:108
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
float ErrElement
Definition: constants.h:55
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
UnitVector3d centroid(VertexIterator const begin, VertexIterator const end)
A base class for image defects.
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
Key< int > photoCalib
Definition: Exposure.cc:67