LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
simpleShape.cc
Go to the documentation of this file.
1// -*- lsst-c++ -*-
2/*
3 * LSST Data Management System
4 * Copyright 2008-2016 AURA/LSST.
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
24#include <array>
25
28#include "lsst/afw/math.h"
29
30namespace lsst { namespace meas { namespace extensions { namespace simpleShape {
31namespace {
32base::FlagDefinitionList flagDefinitions;
33} // end anonymous
34
35base::FlagDefinition const SimpleShape::FAILURE = flagDefinitions.addFailureFlag();
36
38 return flagDefinitions;
39}
40
41
42/* ---- Machinery for iterating over elliptical regions ----------------------------------------------------
43 *
44 * In contrast to FootprintFunctor (which is what's used by most other measurement algorithms), this
45 * machinery doesn't require a virtual function call at each pixel, and hence allows the functor to
46 * be inlined.
47 * Eventually, something like this will replace FootprintFunctor (see #1836)
48 */
49
50namespace {
51
52using afw::geom::Span;
53
54Span clipSpan(Span const & span, geom::Box2I const & box) {
55 if (span.getY() < box.getMinY() || span.getY() > box.getMaxY()) return Span();
56 return Span(span.getY(),
57 std::min(std::max(span.getMinX(), box.getMinX()), box.getMaxX()),
58 std::max(std::min(span.getMaxX(), box.getMaxX()), box.getMinX())
59 );
60}
61
62template <typename Function, typename Iterator>
63void iterateSpan(Function & function, Iterator pixIter, Span const & span) {
64 for (
65 Span::Iterator pointIter = span.begin(), pointEnd = span.end();
66 pointIter != pointEnd;
67 ++pointIter, ++pixIter
68 ) {
69 function(*pointIter, *pixIter);
70 }
71}
72
73template <typename Function, typename Image, typename Region>
74void iterateRegion(Function & function, Image const & image, Region const & region) {
75 geom::Box2I bbox = image.getBBox(afw::image::PARENT);
76 if (bbox.contains(region.getBBox())) {
77 // if the box contains the region, there's no need to check each span to make sure it's entirely
78 // within the image
79 for (
80 typename Region::Iterator spanIter = region.begin(), spanEnd = region.end();
81 spanIter != spanEnd;
82 ++spanIter
83 ) {
84 iterateSpan(
85 function,
86 image.x_at(spanIter->getMinX() - image.getX0(), spanIter->getY() - image.getY0()),
87 *spanIter
88 );
89 }
90 } else {
91 for (
92 typename Region::Iterator spanIter = region.begin(), spanEnd = region.end();
93 spanIter != spanEnd;
94 ++spanIter
95 ) {
96 Span span = clipSpan(*spanIter, bbox);
97 iterateSpan(
98 function,
99 image.x_at(span.getMinX() - image.getX0(), span.getY() - image.getY0()),
100 span
101 );
102 }
103 }
104}
105
106} // anonymous
107
108/* ---- Implementation for SimpleShape algorithm ------------------------------------------------------------
109 *
110 * All of the pixel-based work is done by an accumulating functor that we invoke using the above machinery.
111 * The rest is in the static member function measure(), which provides a way to call the algorithm
112 * outside the context of the pluggable source measurement algorithm framework.
113 */
114
115namespace {
116
117// Enums used to index arrays for code clarity
118enum { Q0=0, QXX, QYY, QXY, QX, QY, N_Q };
119enum { MXX=0, MYY, MXY, MX, MY, N_M };
120
121typedef Eigen::Matrix<double,N_Q,1> VectorQ;
122typedef Eigen::Matrix<double,N_Q,N_Q> MatrixQ;
123typedef Eigen::Matrix<double,N_M,1> VectorM;
124typedef Eigen::Matrix<double,N_M,N_M> MatrixM;
125typedef Eigen::Matrix<double,N_M,N_Q> MatrixMQ;
126
127struct RawMomentAccumulator {
128
129 template <typename PixelT>
130 void operator()(geom::Point2I const & pos, PixelT const & pixel) {
131 geom::Extent2D d = geom::Point2D(pos) - _center;
132 geom::Extent2D gtd = _gt(d);
133 double w = std::exp(-0.5 * (gtd.getX()*gtd.getX() + gtd.getY()*gtd.getY()));
134 VectorQ q = VectorQ::Constant(w);
135 q[QXX] *= d.getX() * d.getX();
136 q[QYY] *= d.getY() * d.getY();
137 q[QXY] *= d.getX() * d.getY();
138 q[QX] *= d.getX();
139 q[QY] *= d.getY();
140 moments += q * pixel.image();
141 covariance.selfadjointView<Eigen::Lower>().rankUpdate(q, pixel.variance());
142 }
143
144 explicit RawMomentAccumulator(afw::geom::ellipses::Ellipse const & weightEllipse) :
145 moments(VectorQ::Zero()),
146 covariance(MatrixQ::Zero()),
147 _center(weightEllipse.getCenter()),
148 _gt(weightEllipse.getCore().getGridTransform())
149 {}
150
151 VectorQ moments;
152 MatrixQ covariance;
153private:
154 geom::Point2D _center;
156};
157
158} // anonymous
159
160template <typename T>
164 double nSigmaRegion
165) {
166 // Define the pixel region we want to accumulate over.
168 regionEllipse.getCore().scale(nSigmaRegion);
169 afw::geom::ellipses::PixelRegion region(regionEllipse);
171
172 // We work in the coordinate system where the origin is the center of the weight function.
173 // This will make things easier in the various transformations we'll have to apply to the raw
174 // moments.
175
176 // All the pixel operations take place in the next two lines, when we use the above machinery
177 // to accumulate the raw moments in a single pass.
178 RawMomentAccumulator functor(weight);
179 iterateRegion(functor, image, region);
180
181 // Then we convert the raw moments to an ellipse and centroid, and propagate the uncertainty
182 // through the covariance matrix.
183 MatrixMQ dm_dq = convertRawMoments(functor.moments, result.ellipse, result.center);
184 MatrixM mCov = dm_dq * functor.covariance.selfadjointView<Eigen::Lower>() * dm_dq.adjoint();
185
186 // Next, we correct the measured moments for the fact that what we just measured was the
187 // moments of the product of the weight function and the data, and we want just the moments
188 // of the data.
189 MatrixM dc_dm = correctWeightedMoments(
191 result.ellipse,
192 result.center
193 );
194 result.covariance = dc_dm * mCov.selfadjointView<Eigen::Lower>() * dc_dm.adjoint();
195
196 // Finally, we switch back to the native image coordinate system.
197 result.center += geom::Extent2D(weight.getCenter());
198
199 return result;
200}
201
202
204 VectorQ const & q,
206 geom::Point2D & center
207) {
208 VectorM m = q.segment<N_M>(1) / q[Q0];
209 MatrixMQ dm_dq = MatrixMQ::Zero();;
210 dm_dq.block<N_M,N_M>(0,1) = MatrixM::Identity() / q[Q0];
211 dm_dq.col(Q0) = -m / q[Q0];
212 // now we account for the centroid^2 terms in the second moments
213 m[MXX] -= m[MX] * m[MX];
214 m[MYY] -= m[MY] * m[MY];
215 m[MXY] -= m[MX] * m[MY];
216 dm_dq(MXX, QX) = -2.0 * m[MX] / q[Q0];
217 dm_dq(MYY, QY) = -2.0 * m[MY] / q[Q0];
218 dm_dq(MXY, QX) = m[MY] / q[Q0];
219 dm_dq(MXY, QY) = m[MX] / q[Q0];
220 double tmp = 2.0 / (q[Q0] * q[Q0] * q[Q0]); // == d(-Q_0^{-2})/d(Q_0)
221 dm_dq(MXX, Q0) += tmp * q[QX] * q[QX];
222 dm_dq(MYY, Q0) += tmp * q[QY] * q[QY];
223 dm_dq(MXY, Q0) += tmp * q[QX] * q[QY];
224
225 ellipse.setIxx(m[MXX]);
226 ellipse.setIyy(m[MYY]);
227 ellipse.setIxy(m[MXY]);
228 center.setX(m[MX]);
229 center.setY(m[MY]);
230
231 return dm_dq;
232}
233
237 geom::Point2D & center
238) {
239 Eigen::Matrix2d wMat = weight.getMatrix();
240 Eigen::Vector2d mVec = center.asEigen();
241 Eigen::Matrix2d mMat = ellipse.getMatrix();
242 if (wMat.determinant() <= 0.0) {
243 throw LSST_EXCEPT(
245 "Weight moments matrix is singular",
247 );
248 }
249 if (mMat.determinant() <= 0.0) {
250 throw LSST_EXCEPT(
252 "Measured moments matrix is singular",
254 );
255 }
256 Eigen::Matrix2d mInv = mMat.inverse();
257 Eigen::Matrix2d cInv = mInv - wMat.inverse();
258 if (cInv.determinant() <= 0.0) {
259 throw LSST_EXCEPT(
261 "Corrected moments matrix is singular",
263 );
264 }
265 Eigen::Matrix2d cMat = cInv.inverse();
266 ellipse.setIxx(cMat(0, 0));
267 ellipse.setIyy(cMat(1, 1));
268 ellipse.setIxy(cMat(0, 1));
269 Eigen::Matrix2d cMat_mInv = cMat * mInv;
270 Eigen::Vector2d mInv_mVec = mInv * mVec;
271 Eigen::Vector2d cVec = cMat_mInv * mVec;
272 center.setX(cVec[0]);
273 center.setY(cVec[1]);
274 Eigen::Matrix2d dcMat_dmxx = cMat_mInv.col(0) * cMat_mInv.col(0).adjoint();
275 Eigen::Matrix2d dcMat_dmyy = cMat_mInv.col(1) * cMat_mInv.col(1).adjoint();
276 Eigen::Matrix2d dcMat_dmxy = cMat_mInv.col(0) * cMat_mInv.col(1).adjoint()
277 + cMat_mInv.col(1) * cMat_mInv.col(0).adjoint();
278 Eigen::Vector2d dcVec_dmxx = cMat_mInv.col(0) * mInv_mVec[0];
279 Eigen::Vector2d dcVec_dmyy = cMat_mInv.col(1) * mInv_mVec[1];
280 Eigen::Vector2d dcVec_dmxy = cMat_mInv.col(0) * mInv_mVec[1] + cMat_mInv.col(1) * mInv_mVec[0];
281 Eigen::Matrix2d const & dcVec_dmVec = cMat_mInv;
282
283 // calculations done - now we just need to put these elements back into a 5x5 matrix to return
284
285 MatrixM dc_dm = MatrixM::Zero();
286 dc_dm(MXX, MXX) = dcMat_dmxx(0, 0);
287 dc_dm(MYY, MXX) = dcMat_dmxx(1, 1);
288 dc_dm(MXY, MXX) = dcMat_dmxx(0, 1);
289 dc_dm(MXX, MYY) = dcMat_dmyy(0, 0);
290 dc_dm(MYY, MYY) = dcMat_dmyy(1, 1);
291 dc_dm(MXY, MYY) = dcMat_dmyy(0, 1);
292 dc_dm(MXX, MXY) = dcMat_dmxy(0, 0);
293 dc_dm(MYY, MXY) = dcMat_dmxy(1, 1);
294 dc_dm(MXY, MXY) = dcMat_dmxy(0, 1);
295 dc_dm(MX, MXX) = dcVec_dmxx[0];
296 dc_dm(MY, MXX) = dcVec_dmxx[1];
297 dc_dm(MX, MYY) = dcVec_dmyy[0];
298 dc_dm(MY, MYY) = dcVec_dmyy[1];
299 dc_dm(MX, MXY) = dcVec_dmxy[0];
300 dc_dm(MY, MXY) = dcVec_dmxy[1];
301 dc_dm(MX, MX) = dcVec_dmVec(0, 0);
302 dc_dm(MX, MY) = dcVec_dmVec(0, 1);
303 dc_dm(MY, MX) = dcVec_dmVec(1, 0);
304 dc_dm(MY, MY) = dcVec_dmVec(1, 1);
305
306 return dc_dm;
307}
308
309
310SimpleShapeResult::SimpleShapeResult() : ellipse(std::numeric_limits<lsst::meas::base::ErrElement>::quiet_NaN(),
311 std::numeric_limits<lsst::meas::base::ErrElement>::quiet_NaN(),
312 std::numeric_limits<lsst::meas::base::ErrElement>::quiet_NaN()),
313 center(std::numeric_limits<lsst::meas::base::ErrElement>::quiet_NaN(),
314 std::numeric_limits<lsst::meas::base::ErrElement>::quiet_NaN()),
315 covariance(Eigen::Matrix<double,5,5>::Constant(std::numeric_limits<lsst::meas::base::ErrElement>::quiet_NaN()))
316{}
317
318
320 afw::table::Schema & schema,
321 std::string const & name
322) {
324 r._shapeResult = lsst::afw::table::QuadrupoleKey::addFields(schema,
325 name, "elliptical Gaussian moments");
326 r._centroidResult = lsst::afw::table::Point2DKey::addFields(schema,
327 name, "elliptical Gaussian moments", "pixel");
328 r._uncertantyResult = lsst::afw::table::CovarianceMatrixKey<double, 5>::addFields(schema,name,
329 std::vector<std::string> ({"Ixx", "Iyy", "Ixy",
330 "Ix", "Iy"}), "pixel");
331 r._flagHandler = lsst::meas::base::FlagHandler::addFields(schema,
333 return r;
334}
335
337 _shapeResult(s),
338 _centroidResult(s),
339 _uncertantyResult(s, std::vector<std::string> ({"Ixx", "Iyy", "Ixy", "Ix", "Iy"})),
340 _flagHandler(s, SimpleShape::getFlagDefinitions())
341{}
342
345 result.ellipse = record.get(_shapeResult);
346 result.center = record.get(_centroidResult);
347 result.covariance = record.get(_uncertantyResult);
348 for (std::size_t n = 0; n < SimpleShape::N_FLAGS; ++n) {
349 result.flags[n] = _flagHandler.getValue(record, n);
350 }
351 return result;
352}
353
355 record.set(_shapeResult, value.ellipse);
356 record.set(_centroidResult, value.center);
357 record.set(_uncertantyResult, value.covariance);
358 for (std::size_t n = 0; n < SimpleShape::N_FLAGS; ++n) {
359 _flagHandler.setValue(record, n, value.flags[n]);
360 }
361}
362
364 return _shapeResult == other._shapeResult &&
365 _centroidResult == other._centroidResult &&
366 _uncertantyResult == other._uncertantyResult;
367 //don't bother with flags - if we've gotten this far, it's basically impossible the flags don't match
368}
369
371 return _shapeResult.isValid() &&
372 _centroidResult.isValid() &&
373 _uncertantyResult.isValid();
374 //don't bother with flags - if we've gotten this far, it's basically impossible the flags don't match
375}
376
377
379 Control const & ctrl,
380 std::string const & name,
381 afw::table::Schema & schema
382 )
383 : _ctrl(ctrl),
384 _resultKey(SimpleShapeResultKey::addFields(schema, name)),
385 _centroidExtractor(schema, name)
386 {}
387
390 afw::image::Exposure<float> const & exposure
391) const {
392 geom::Point2D center = _centroidExtractor(source, _resultKey.getFlagHandler());
394 // set flags so an exception throw produces a flagged source
395 SimpleShapeResult result = computeMoments(weight, exposure.getMaskedImage(), _ctrl.nSigmaRegion);
396 source.set(_resultKey, result);
397}
398
400 afw::table::SourceRecord & measRecord,
402) const {
403 _resultKey.getFlagHandler().handleFailure(measRecord, error);
404}
405
406#define INSTANTIATE_IMAGE(IMAGE) \
407 template SimpleShapeResult SimpleShape::computeMoments(\
408 afw::geom::ellipses::Ellipse const &, \
409 IMAGE const &, \
410 double \
411 )
412
416
417}}}} // namespace lsst::meas::extensions::simpleShape
py::object result
Definition _schema.cc:429
AmpInfoBoxKey bbox
Definition Amplifier.cc:117
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
int m
Definition SpanSet.cc:48
A range of pixels within one row of an Image.
Definition Span.h:47
An ellipse core for the semimajor/semiminor axis and position angle parametrization (a,...
Definition Axes.h:47
void scale(double factor)
Scale the size of the ellipse core by the given factor.
Definition BaseCore.cc:103
An ellipse defined by an arbitrary BaseCore and a center point.
Definition Ellipse.h:51
BaseCore const & getCore() const
Return the ellipse core.
Definition Ellipse.h:71
A pixelized region containing all pixels whose centers are within an Ellipse.
Definition PixelRegion.h:46
An ellipse core with quadrupole moments as parameters.
Definition Quadrupole.h:47
Matrix const & getMatrix() const
Return a 2x2 symmetric matrix of the parameters.
Definition Quadrupole.h:80
A class to contain the data, WCS, and other information needed to describe an image of the sky.
Definition Exposure.h:72
A class to manipulate images, masks, and variance as a single object.
Definition MaskedImage.h:74
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
static CovarianceMatrixKey addFields(Schema &schema, std::string const &prefix, NameArray const &names, std::string const &unit, bool diagonalOnly=false)
Add covariance matrix fields to a Schema, and return a CovarianceMatrixKey to manage 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
bool isValid() const noexcept
Return True if both the x and y Keys are valid.
Definition aggregates.h:106
static QuadrupoleKey addFields(Schema &schema, std::string const &name, std::string const &doc, CoordinateType coordType=CoordinateType::PIXEL)
Add a set of quadrupole subfields to a schema and return a QuadrupoleKey that points to them.
bool isValid() const noexcept
Return True if all the constituent Keys are valid.
Definition aggregates.h:429
Defines the fields and offsets for a table.
Definition Schema.h:51
Record class that contains measurements made on a single exposure.
Definition Source.h:78
A proxy type for name lookups in a Schema.
Definition Schema.h:367
An integer coordinate rectangle.
Definition Box.h:55
int getMinY() const noexcept
Definition Box.h:158
int getMinX() const noexcept
Definition Box.h:157
int getMaxX() const noexcept
Definition Box.h:161
int getMaxY() const noexcept
Definition Box.h:162
EigenVector const & asEigen() const noexcept(IS_ELEMENT_NOTHROW_COPYABLE)
Return a fixed-size Eigen representation of the coordinate object.
A 2D linear coordinate transformation.
vector-type utility class to build a collection of FlagDefinitions
Definition FlagHandler.h:60
void handleFailure(afw::table::BaseRecord &record, MeasurementError const *error=nullptr) const
Handle an expected or unexpected Exception thrown by a measurement algorithm.
static FlagHandler addFields(afw::table::Schema &schema, std::string const &prefix, FlagDefinitionList const &flagDefs, FlagDefinitionList const &exclDefs=FlagDefinitionList::getEmptyList())
Add Flag fields to a schema, creating a FlagHandler object to manage them.
void setValue(afw::table::BaseRecord &record, std::size_t i, bool value) const
Set the flag field corresponding to the given flag index.
bool getValue(afw::table::BaseRecord const &record, std::size_t i) const
Return the value of the flag field corresponding to the given flag index.
Exception to be thrown when a measurement algorithm experiences a known failure mode.
Definition exceptions.h:48
A C++ control class to handle SdssShapeAlgorithm's configuration.
Definition simpleShape.h:44
double sigma
"Sigma of circular Gaussian used as weight function, in pixels" ;
Definition simpleShape.h:46
double nSigmaRegion
"Maximum radius for pixels to include, in units of sigma" ;
Definition simpleShape.h:47
SimpleShape(Control const &ctrl, std::string const &name, afw::table::Schema &schema)
static Eigen::Matrix< double, 5, 5 > correctWeightedMoments(afw::geom::ellipses::Quadrupole const &weight, afw::geom::ellipses::Quadrupole &ellipse, geom::Point2D &center)
Correct moments measured with a Gaussian weight function by assuming the data was also an elliptical ...
static SimpleShapeResult computeMoments(afw::geom::ellipses::Ellipse const &weight, afw::image::MaskedImage< T > const &image, double nSigmaRegion=3.0)
Compute the Gaussian-weighted moments of an image.
virtual void fail(afw::table::SourceRecord &measRecord, lsst::meas::base::MeasurementError *error=NULL) const
Handle an exception thrown by the current algorithm by setting flags in the given record.
static base::FlagDefinitionList const & getFlagDefinitions()
static base::FlagDefinition const FAILURE
static Eigen::Matrix< double, 5, 6 > convertRawMoments(Eigen::Matrix< double, 6, 1 > const &q, afw::geom::ellipses::Quadrupole &quadrupole, geom::Point2D &center)
Convert linear raw moments into an ellipse and centroid, and return the derivative of the conversion.
virtual void measure(afw::table::SourceRecord &measRecord, afw::image::Exposure< float > const &exposure) const
Called to measure a single child source in an image.
Struct to hold the results of SimpleShape when we don't run it as a plugin.
SimpleShapeResult()
Constructor; initializes everything to Nan.
afw::geom::ellipses::Quadrupole ellipse
Measured second moments.
virtual void set(afw::table::BaseRecord &record, SimpleShapeResult const &value) const
virtual SimpleShapeResult get(afw::table::BaseRecord const &record) const
bool isValid() const
Return True if the key is valid.
bool operator==(SimpleShapeResultKey const &other) const
Compare the FunctorKey for equality with another, using the underlying Keys.
lsst::meas::base::FlagHandler const & getFlagHandler() const
Definition simpleShape.h:85
static SimpleShapeResultKey addFields(afw::table::Schema &schema, std::string const &name)
T exp(T... args)
T max(T... args)
#define INSTANTIATE_IMAGE(IMAGE)
Definition ImagePca.cc:125
T min(T... args)
Extent< double, 2 > Extent2D
Definition Extent.h:400
Point< double, 2 > Point2D
Definition Point.h:324
STL namespace.
afw::table::Key< double > weight
double w
Definition CoaddPsf.cc:69
MatrixQ covariance
VectorQ moments