lsst::meas::extensions::simpleShape::SimpleShape Class Reference

#include <simpleShape.h>

Inheritance diagram for lsst::meas::extensions::simpleShape::SimpleShape:

Public Types
typedef SimpleShapeControl	Control

Public Member Functions
	SimpleShape (Control const &ctrl, std::string const &name, afw::table::Schema &schema)
virtual void	measure (afw::table::SourceRecord &measRecord, afw::image::Exposure< float > const &exposure) const
	Called to measure a single child source in an image. More...
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. More...
virtual void	measureForced (afw::table::SourceRecord &measRecord, afw::image::Exposure< float > const &exposure, afw::table::SourceRecord const &refRecord, afw::geom::SkyWcs const &refWcs) const
	Called to measure a single child source in an image. More...
virtual void	measureNForced (afw::table::SourceCatalog const &measCat, afw::image::Exposure< float > const &exposure, afw::table::SourceCatalog const &refRecord, afw::geom::SkyWcs const &refWcs) const
	Called to simultaneously measure all children in a deblend family, in a single image. More...
virtual void	measureN (afw::table::SourceCatalog const &measCat, afw::image::Exposure< float > const &exposure) const
	Called to simultaneously measure all children in a deblend family, in a single image. More...
std::string	getLogName () const

Static Public Member Functions
static base::FlagDefinitionList const &	getFlagDefinitions ()
template<typename T >
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. More...
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. More...
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 Gaussian, and return the derivative of the correction. More...

Static Public Attributes
static unsigned int const	N_FLAGS = 1
static base::FlagDefinition const	FAILURE = flagDefinitions.addFailureFlag()

Protected Attributes
std::string	_logName

Detailed Description

Definition at line 95 of file simpleShape.h.

Member Typedef Documentation

Control

typedef SimpleShapeControl lsst::meas::extensions::simpleShape::SimpleShape::Control

Definition at line 103 of file simpleShape.h.

Constructor & Destructor Documentation

SimpleShape()

lsst::meas::extensions::simpleShape::SimpleShape::SimpleShape	(	Control const &	ctrl,
		std::string const &	name,
		afw::table::Schema &	schema
	)

Definition at line 378 of file simpleShape.cc.

    : _ctrl(ctrl),
      _resultKey(SimpleShapeResultKey::addFields(schema, name)),
      _centroidExtractor(schema, name)
    {}

Member Function Documentation

computeMoments()

template<typename T >

template SimpleShapeResult lsst::meas::extensions::simpleShape::SimpleShape::computeMoments ( afw::geom::ellipses::Ellipse const &  weight, afw::image::MaskedImage< T > const &  image, double  nSigmaRegion = 3.0  )

static

Compute the Gaussian-weighted moments of an image.

Parameters

[in]	weight	An ellipse object of Gaussian weights to apply to the measurement.
[in]	image	A Masked image instance with int float or double pixels.
[in]	nSigmaRegion	Maximum radius for pixels to include, in units of sigma

Definition at line 161 of file simpleShape.cc.

  {
    // Define the pixel region we want to accumulate over.
    afw::geom::ellipses::Ellipse regionEllipse(weight);
    regionEllipse.getCore().scale(nSigmaRegion);
    afw::geom::ellipses::PixelRegion region(regionEllipse);
    SimpleShapeResult result;
 
    // We work in the coordinate system where the origin is the center of the weight function.
    // This will make things easier in the various transformations we'll have to apply to the raw
    // moments.
 
    // All the pixel operations take place in the next two lines, when we use the above machinery
    // to accumulate the raw moments in a single pass.
    RawMomentAccumulator functor(weight);
    iterateRegion(functor, image, region);
 
    // Then we convert the raw moments to an ellipse and centroid, and propagate the uncertainty
    // through the covariance matrix.
    MatrixMQ dm_dq = convertRawMoments(functor.moments, result.ellipse, result.center);
    MatrixM mCov = dm_dq * functor.covariance.selfadjointView<Eigen::Lower>() * dm_dq.adjoint();
 
    // Next, we correct the measured moments for the fact that what we just measured was the
    // moments of the product of the weight function and the data, and we want just the moments
    // of the data.
    MatrixM dc_dm = correctWeightedMoments(
        weight.getCore().as<afw::geom::ellipses::Quadrupole>(),
        result.ellipse,
        result.center
    );
    result.covariance = dc_dm * mCov.selfadjointView<Eigen::Lower>() * dc_dm.adjoint();
 
    // Finally, we switch back to the native image coordinate system.
    result.center += geom::Extent2D(weight.getCenter());
 
    return result;
}

convertRawMoments()

MatrixMQ lsst::meas::extensions::simpleShape::SimpleShape::convertRawMoments ( Eigen::Matrix< double, 6, 1 > const &  q, afw::geom::ellipses::Quadrupole &  quadrupole, geom::Point2D &  center  )

static

Convert linear raw moments into an ellipse and centroid, and return the derivative of the conversion.

Note

This function is mainly intended for internal use, and is only exposed publically so it can be unit-tested in Python.

For weight function \(w\) and data \(p\), the "raw" moments \(Q\) are defined as:

\begin{eqnarray*} Q_0 &=& \sum_n w(x_n, y_n) p_n \\ Q_{xx} &=& \sum_n w(x_n, y_n) x_n^2 p_n \\ Q_{yy} &=& \sum_n w(x_n, y_n) y_n^2 p_n \\ Q_{xy} &=& \sum_n w(x_n, y_n) x_n y_n p_n \\ Q_x &=& \sum_n w(x_n, y_n) x_n p_n \\ Q_y &=& \sum_n w(x_n, y_n) y_n p_n \end{eqnarray*}

whereas the converted ellipse and centroid moments are:

\begin{eqnarray*} M_{xx} &=& Q_{xx} / Q_0 - Q_x^2 \\ M_{xx} &=& Q_{yy} / Q_0 - Q_y^2 \\ M_{xx} &=& Q_{xy} / Q_0 - Q_x Q_y \\ M_x &=& Q_x / Q_0 \\ M_y &=& Q_y / Q_0 \end{eqnarray*}

Note the slightly unusual ordering; this is for consistency with afw::geom::ellipses::Ellipse.

Definition at line 203 of file simpleShape.cc.

  {
    VectorM m = q.segment<N_M>(1) / q[Q0];
    MatrixMQ dm_dq = MatrixMQ::Zero();;
    dm_dq.block<N_M,N_M>(0,1) = MatrixM::Identity() / q[Q0];
    dm_dq.col(Q0) = -m / q[Q0];
    // now we account for the centroid^2 terms in the second moments
    m[MXX] -= m[MX] * m[MX];
    m[MYY] -= m[MY] * m[MY];
    m[MXY] -= m[MX] * m[MY];
    dm_dq(MXX, QX) = -2.0 * m[MX] / q[Q0];
    dm_dq(MYY, QY) = -2.0 * m[MY] / q[Q0];
    dm_dq(MXY, QX) = m[MY] / q[Q0];
    dm_dq(MXY, QY) = m[MX] / q[Q0];
    double tmp = 2.0 / (q[Q0] * q[Q0] * q[Q0]); // == d(-Q_0^{-2})/d(Q_0)
    dm_dq(MXX, Q0) += tmp * q[QX] * q[QX];
    dm_dq(MYY, Q0) += tmp * q[QY] * q[QY];
    dm_dq(MXY, Q0) += tmp * q[QX] * q[QY];
 
    ellipse.setIxx(m[MXX]);
    ellipse.setIyy(m[MYY]);
    ellipse.setIxy(m[MXY]);
    center.setX(m[MX]);
    center.setY(m[MY]);
 
    return dm_dq;
}

correctWeightedMoments()

MatrixM lsst::meas::extensions::simpleShape::SimpleShape::correctWeightedMoments ( afw::geom::ellipses::Quadrupole const &  weight, afw::geom::ellipses::Quadrupole &  ellipse, geom::Point2D &  center  )

static

Correct moments measured with a Gaussian weight function by assuming the data was also an elliptical Gaussian, and return the derivative of the correction.

Note

This function is mainly intended for internal use, and is only exposed publically so it can be unit-tested in Python.

If we naively measure Gaussian-weighted moments, we'll measure the moments of the product of the weight function and the data. What we want is the moments of the data, as if we had measured them with no weight function (but without sacrificing the S/N benefit that comes from using a weight function). To do that, we assume the data is also an elliptical Gaussian, and "divide" the weight function from the measured moments to compute it.

If \(W\) and \(M\) are the quadruple matrices of the weight function and measurement, and \(\eta\) is the measured centroid (we work in a coordinate system where the weight function is centered at the origin), then the corrected quadrupole matrix \(C\) and centroid are \(\nu\) are:

\begin{eqnarray*} C &=& \left(M^{-1} - W^{-1}\right)^{-1} \\ \nu &=& C M^{-1} \eta \end{eqnarray*}

Definition at line 234 of file simpleShape.cc.

  {
    Eigen::Matrix2d wMat = weight.getMatrix();
    Eigen::Vector2d mVec = center.asEigen();
    Eigen::Matrix2d mMat = ellipse.getMatrix();
    if (wMat.determinant() <= 0.0) {
        throw LSST_EXCEPT(
            base::MeasurementError,
            "Weight moments matrix is singular",
            FAILURE.number
        );
    }
    if (mMat.determinant() <= 0.0) {
        throw LSST_EXCEPT(
            base::MeasurementError,
            "Measured moments matrix is singular",
            FAILURE.number
        );
    }
    Eigen::Matrix2d mInv = mMat.inverse();
    Eigen::Matrix2d cInv = mInv - wMat.inverse();
    if (cInv.determinant() <= 0.0) {
        throw LSST_EXCEPT(
            base::MeasurementError,
            "Corrected moments matrix is singular",
            FAILURE.number
        );
    }
    Eigen::Matrix2d cMat = cInv.inverse();
    ellipse.setIxx(cMat(0, 0));
    ellipse.setIyy(cMat(1, 1));
    ellipse.setIxy(cMat(0, 1));
    Eigen::Matrix2d cMat_mInv = cMat * mInv;
    Eigen::Vector2d mInv_mVec = mInv * mVec;
    Eigen::Vector2d cVec = cMat_mInv * mVec;
    center.setX(cVec[0]);
    center.setY(cVec[1]);
    Eigen::Matrix2d dcMat_dmxx = cMat_mInv.col(0) * cMat_mInv.col(0).adjoint();
    Eigen::Matrix2d dcMat_dmyy = cMat_mInv.col(1) * cMat_mInv.col(1).adjoint();
    Eigen::Matrix2d dcMat_dmxy = cMat_mInv.col(0) * cMat_mInv.col(1).adjoint()
        + cMat_mInv.col(1) * cMat_mInv.col(0).adjoint();
    Eigen::Vector2d dcVec_dmxx = cMat_mInv.col(0) * mInv_mVec[0];
    Eigen::Vector2d dcVec_dmyy = cMat_mInv.col(1) * mInv_mVec[1];
    Eigen::Vector2d dcVec_dmxy = cMat_mInv.col(0) * mInv_mVec[1] + cMat_mInv.col(1) * mInv_mVec[0];
    Eigen::Matrix2d const & dcVec_dmVec = cMat_mInv;
 
    // calculations done - now we just need to put these elements back into a 5x5 matrix to return
 
    MatrixM dc_dm = MatrixM::Zero();
    dc_dm(MXX, MXX) = dcMat_dmxx(0, 0);
    dc_dm(MYY, MXX) = dcMat_dmxx(1, 1);
    dc_dm(MXY, MXX) = dcMat_dmxx(0, 1);
    dc_dm(MXX, MYY) = dcMat_dmyy(0, 0);
    dc_dm(MYY, MYY) = dcMat_dmyy(1, 1);
    dc_dm(MXY, MYY) = dcMat_dmyy(0, 1);
    dc_dm(MXX, MXY) = dcMat_dmxy(0, 0);
    dc_dm(MYY, MXY) = dcMat_dmxy(1, 1);
    dc_dm(MXY, MXY) = dcMat_dmxy(0, 1);
    dc_dm(MX, MXX) = dcVec_dmxx[0];
    dc_dm(MY, MXX) = dcVec_dmxx[1];
    dc_dm(MX, MYY) = dcVec_dmyy[0];
    dc_dm(MY, MYY) = dcVec_dmyy[1];
    dc_dm(MX, MXY) = dcVec_dmxy[0];
    dc_dm(MY, MXY) = dcVec_dmxy[1];
    dc_dm(MX, MX) = dcVec_dmVec(0, 0);
    dc_dm(MX, MY) = dcVec_dmVec(0, 1);
    dc_dm(MY, MX) = dcVec_dmVec(1, 0);
    dc_dm(MY, MY) = dcVec_dmVec(1, 1);
 
    return dc_dm;
}

fail()

void lsst::meas::extensions::simpleShape::SimpleShape::fail ( afw::table::SourceRecord &  measRecord, lsst::meas::base::MeasurementError *  error = NULL  ) const

virtual

Handle an exception thrown by the current algorithm by setting flags in the given record.

fail() is called by the measurement framework when an exception is allowed to propagate out of one the algorithm's measure() methods. It should generally set both a general failure flag for the algorithm as well as a specific flag indicating the error condition, if possible. To aid in this, if the exception was an instance of MeasurementError, it will be passed in, carrying information about what flag to set.

An algorithm can also to chose to set flags within its own measure() methods, and then just return, rather than throw an exception. However, fail() should be implemented even when all known failure modes do not throw exceptions, to ensure that unexpected exceptions thrown in lower-level code are properly handled.

Implements lsst::meas::base::BaseAlgorithm.

Definition at line 399 of file simpleShape.cc.

        {
    _resultKey.getFlagHandler().handleFailure(measRecord, error);
}

getFlagDefinitions()

base::FlagDefinitionList const & lsst::meas::extensions::simpleShape::SimpleShape::getFlagDefinitions ( )

static

Definition at line 37 of file simpleShape.cc.

                                                               {
    return flagDefinitions;
}

getLogName()

std::string lsst::meas::base::BaseAlgorithm::getLogName ( ) const

inlineinherited

Definition at line 66 of file Algorithm.h.

{ return _logName; }

measure()

void lsst::meas::extensions::simpleShape::SimpleShape::measure ( afw::table::SourceRecord &  measRecord, afw::image::Exposure< float > const &  exposure  ) const

virtual

Called to measure a single child source in an image.

Before this method is called, all neighbors will be replaced with noise, using the outputs of the deblender. Outputs should be saved in the given SourceRecord, which can also be used to obtain centroid (see SafeCentroidExtractor) and shape (see SafeShapeExtractor) information.

Implements lsst::meas::base::SingleFrameAlgorithm.

Definition at line 388 of file simpleShape.cc.

        {
    geom::Point2D center = _centroidExtractor(source, _resultKey.getFlagHandler());
    afw::geom::ellipses::Ellipse weight(afw::geom::ellipses::Axes(_ctrl.sigma), center);
    // set flags so an exception throw produces a flagged source
    SimpleShapeResult result = computeMoments(weight, exposure.getMaskedImage(), _ctrl.nSigmaRegion);
    source.set(_resultKey, result);
}

measureForced()

virtual void lsst::meas::base::SimpleAlgorithm::measureForced ( afw::table::SourceRecord &  measRecord, afw::image::Exposure< float > const &  exposure, afw::table::SourceRecord const &  refRecord, afw::geom::SkyWcs const &  refWcs  ) const

inlinevirtualinherited

Called to measure a single child source in an image.

Before this method is called, all neighbors will be replaced with noise, using the outputs of the deblender. Outputs should be saved in the given SourceRecord, which can also be used to obtain centroid (see SafeCentroidExtractor) and shape (see SafeShapeExtractor) information.

Implements lsst::meas::base::ForcedAlgorithm.

Reimplemented in lsst::meas::extensions::photometryKron::KronFluxAlgorithm.

Definition at line 172 of file Algorithm.h.

                                                                  {
        measure(measRecord, exposure);
    }

measureN()

void lsst::meas::base::SingleFrameAlgorithm::measureN ( afw::table::SourceCatalog const &  measCat, afw::image::Exposure< float > const &  exposure  ) const

virtualinherited

Called to simultaneously measure all children in a deblend family, in a single image.

Outputs should be saved in the given SourceCatalog, which can also be used to obtain centroid (see SafeCentroidExtractor) and shape (see SafeShapeExtractor) information.

The default implementation simply throws an exception, indicating that simultaneous measurement is not supported.

Definition at line 31 of file Algorithm.cc.

                                                                                   {
    throw LSST_EXCEPT(pex::exceptions::LogicError, "measureN not implemented for this algorithm");
}

measureNForced()

virtual void lsst::meas::base::SimpleAlgorithm::measureNForced ( afw::table::SourceCatalog const &  measCat, afw::image::Exposure< float > const &  exposure, afw::table::SourceCatalog const &  refRecord, afw::geom::SkyWcs const &  refWcs  ) const

inlinevirtualinherited

Called to simultaneously measure all children in a deblend family, in a single image.

Outputs should be saved in the given SourceCatalog, which can also be used to obtain centroid (see SafeCentroidExtractor) and shape (see SafeShapeExtractor) information.

The default implementation simply throws an exception, indicating that simultaneous measurement is not supported.

Reimplemented from lsst::meas::base::ForcedAlgorithm.

Definition at line 179 of file Algorithm.h.

                                                                   {
        measureN(measCat, exposure);
    }

Member Data Documentation

_logName

std::string lsst::meas::base::BaseAlgorithm::_logName

protectedinherited

Definition at line 69 of file Algorithm.h.

FAILURE

base::FlagDefinition const lsst::meas::extensions::simpleShape::SimpleShape::FAILURE = flagDefinitions.addFailureFlag()

static

Definition at line 101 of file simpleShape.h.

N_FLAGS

unsigned int const lsst::meas::extensions::simpleShape::SimpleShape::N_FLAGS = 1

static

Definition at line 100 of file simpleShape.h.

---

The documentation for this class was generated from the following files:

• /j/snowflake/release/lsstsw/stack/lsst-scipipe-0.7.0/Linux64/meas_extensions_simpleShape/22.0.0+8d77f4f51a/include/lsst/meas/extensions/simpleShape.h
• /j/snowflake/release/lsstsw/stack/lsst-scipipe-0.7.0/Linux64/meas_extensions_simpleShape/22.0.0+8d77f4f51a/src/simpleShape.cc