LSSTApplications
20.0.0
LSSTDataManagementBasePackage
|
Go to the documentation of this file.
30 namespace lsst {
namespace meas {
namespace extensions {
namespace simpleShape {
32 base::FlagDefinitionList flagDefinitions;
38 return flagDefinitions;
54 Span 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(),
62 template <
typename Function,
typename Iterator>
63 void iterateSpan(Function &
function,
Iterator pixIter, Span
const & span) {
66 pointIter != pointEnd;
67 ++pointIter, ++pixIter
69 function(*pointIter, *pixIter);
73 template <
typename Function,
typename Image,
typename Region>
74 void iterateRegion(Function &
function, Image
const &
image, Region
const & region) {
76 if (
bbox.contains(region.getBBox())) {
86 image.x_at(spanIter->getMinX() -
image.getX0(), spanIter->getY() -
image.getY0()),
96 Span span = clipSpan(*spanIter,
bbox);
99 image.x_at(span.getMinX() -
image.getX0(), span.getY() -
image.getY0()),
118 enum { Q0=0, QXX, QYY, QXY, QX, QY, N_Q };
119 enum { MXX=0, MYY, MXY, MX, MY, N_M };
121 typedef Eigen::Matrix<double,N_Q,1> VectorQ;
122 typedef Eigen::Matrix<double,N_Q,N_Q> MatrixQ;
123 typedef Eigen::Matrix<double,N_M,1> VectorM;
124 typedef Eigen::Matrix<double,N_M,N_M> MatrixM;
125 typedef Eigen::Matrix<double,N_M,N_Q> MatrixMQ;
127 struct RawMomentAccumulator {
129 template <
typename PixelT>
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();
141 covariance.selfadjointView<Eigen::Lower>().rankUpdate(q,
pixel.variance());
144 explicit RawMomentAccumulator(afw::geom::ellipses::Ellipse
const & weightEllipse) :
147 _center(weightEllipse.getCenter()),
148 _gt(weightEllipse.getCore().getGridTransform())
160 template <
typename T>
178 RawMomentAccumulator functor(
weight);
179 iterateRegion(functor,
image, region);
184 MatrixM mCov = dm_dq * functor.covariance.selfadjointView<Eigen::Lower>() * dm_dq.adjoint();
194 result.covariance = dc_dm * mCov.selfadjointView<Eigen::Lower>() * dc_dm.adjoint();
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];
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]);
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];
239 Eigen::Matrix2d wMat =
weight.getMatrix();
240 Eigen::Vector2d mVec = center.asEigen();
241 Eigen::Matrix2d mMat = ellipse.
getMatrix();
242 if (wMat.determinant() <= 0.0) {
245 "Weight moments matrix is singular",
249 if (mMat.determinant() <= 0.0) {
252 "Measured moments matrix is singular",
256 Eigen::Matrix2d mInv = mMat.inverse();
257 Eigen::Matrix2d cInv = mInv - wMat.inverse();
258 if (cInv.determinant() <= 0.0) {
261 "Corrected moments matrix is singular",
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;
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);
325 name,
"elliptical Gaussian moments");
327 name,
"elliptical Gaussian moments",
"pixel");
330 "Ix",
"Iy"}),
"pixel");
339 _uncertantyResult(s,
std::vector<
std::string> ({
"Ixx",
"Iyy",
"Ixy",
"Ix",
"Iy"})),
345 result.ellipse = record.
get(_shapeResult);
346 result.center = record.
get(_centroidResult);
347 result.covariance = record.
get(_uncertantyResult);
356 record.
set(_centroidResult, value.
center);
364 return _shapeResult ==
other._shapeResult &&
365 _centroidResult ==
other._centroidResult &&
366 _uncertantyResult ==
other._uncertantyResult;
371 return _shapeResult.
isValid() &&
406 #define INSTANTIATE_IMAGE(IMAGE) \
407 template SimpleShapeResult SimpleShape::computeMoments(\
408 afw::geom::ellipses::Ellipse const &, \
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
A range of pixels within one row of an Image.
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 ¢er)
Convert linear raw moments into an ellipse and centroid, and return the derivative of the conversion.
void handleFailure(afw::table::BaseRecord &record, MeasurementError const *error=nullptr) const
Handle an expected or unexpected Exception thrown by a measurement algorithm.
Record class that contains measurements made on a single exposure.
Field< T >::Value get(Key< T > const &key) const
Return the value of a field for the given key.
double nSigmaRegion
"Maximum radius for pixels to include, in units of sigma" ;
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.
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.
virtual SimpleShapeResult get(afw::table::BaseRecord const &record) const
bool isValid() const
Return True if the key is valid.
A C++ control class to handle SdssShapeAlgorithm's configuration.
Struct to hold the results of SimpleShape when we don't run it as a plugin.
lsst::meas::base::FlagHandler const & getFlagHandler() const
Exception to be thrown when a measurement algorithm experiences a known failure mode.
std::bitset< SimpleShape::N_FLAGS > flags
Defines the fields and offsets for a table.
An ellipse core for the semimajor/semiminor axis and position angle parametrization (a,...
Point< double, 2 > Point2D
Eigen::Matrix< double, 5, 5 > covariance
Matrix of uncertainties; ordered Ixx, Iyy, Ixy, Ix, Iy.
MaskedImageT getMaskedImage()
Return the MaskedImage.
Matrix const & getMatrix() const
Return a 2x2 symmetric matrix of the parameters.
vector-type utility class to build a collection of FlagDefinitions
afw::geom::ellipses::Quadrupole ellipse
Measured second moments.
bool isValid() const noexcept
Return True if all the constituent error Keys are valid.
bool isValid() const noexcept
Return True if both the x and y Keys are valid.
static Eigen::Matrix< double, 5, 5 > correctWeightedMoments(afw::geom::ellipses::Quadrupole const &weight, afw::geom::ellipses::Quadrupole &ellipse, geom::Point2D ¢er)
Correct moments measured with a Gaussian weight function by assuming the data was also an elliptical ...
void setValue(afw::table::BaseRecord &record, std::size_t i, bool value) const
Set the flag field corresponding to the given flag index.
static unsigned int const N_FLAGS
A class to manipulate images, masks, and variance as a single object.
const char * source()
Source function that allows astChannel to source from a Stream.
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.
ItemVariant const * other
Base class for all records.
geom::Point2D center
Measured first moments, or the input center if !recentroid.
FastFinder::Iterator Iterator
SimpleShape(Control const &ctrl, std::string const &name, afw::table::Schema &schema)
virtual void set(afw::table::BaseRecord &record, SimpleShapeResult const &value) const
table::PointKey< int > pixel
afw::table::Key< double > weight
A base class for image defects.
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
An ellipse defined by an arbitrary BaseCore and a center point.
A pixelized region containing all pixels whose centers are within an Ellipse.
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > Matrix
An ellipse core with quadrupole moments as parameters.
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.
int getMaxY() const noexcept
A proxy type for name lookups in a Schema.
int getMaxX() const noexcept
An integer coordinate rectangle.
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.
int getMinX() const noexcept
Extent< double, 2 > Extent2D
bool operator==(SimpleShapeResultKey const &other) const
Compare the FunctorKey for equality with another, using the underlying Keys.
virtual void measure(afw::table::SourceRecord &measRecord, afw::image::Exposure< float > const &exposure) const
Called to measure a single child source in 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.
BaseCore const & getCore() const
Return the ellipse core.
void set(Key< T > const &key, U const &value)
Set value of a field for the given key.
SimpleShapeResult()
Constructor; initializes everything to Nan.
static base::FlagDefinitionList const & getFlagDefinitions()
void scale(double factor)
Scale the size of the ellipse core by the given factor.
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.
double sigma
"Sigma of circular Gaussian used as weight function, in pixels" ;
int getMinY() const noexcept
static SimpleShapeResultKey addFields(afw::table::Schema &schema, std::string const &name)
bool isValid() const noexcept
Return True if all the constituent Keys are valid.
#define INSTANTIATE_IMAGE(IMAGE)