lsst::ip::diffim::PsfDipoleFlux Class Reference

#include <DipoleAlgorithms.h>

Inheritance diagram for lsst::ip::diffim::PsfDipoleFlux:

Public Types
typedef PsfDipoleFluxControl	Control
Public Types inherited from lsst::ip::diffim::DipoleFluxAlgorithm
enum	{ FAILURE =meas::base::FlagHandler::FAILURE, POS_FAILURE, NEG_FAILURE, N_FLAGS }
typedef DipoleFluxControl	Control
typedef meas::base::FluxResultKey	ResultKey

Public Member Functions
	PsfDipoleFlux (PsfDipoleFluxControl const &ctrl, std::string const &name, afw::table::Schema &schema)
std::pair< double, int >	chi2 (afw::table::SourceRecord &source, afw::image::Exposure< float > const &exposure, double negCenterX, double negCenterY, double negFlux, double posCenterX, double poCenterY, double posFlux) const
void	measure (afw::table::SourceRecord &measRecord, afw::image::Exposure< float > const &exposure) const
void	fail (afw::table::SourceRecord &measRecord, meas::base::MeasurementError *error=NULL) const
Public Member Functions inherited from lsst::ip::diffim::DipoleFluxAlgorithm
	DipoleFluxAlgorithm (Control const &ctrl, std::string const &name, afw::table::Schema &schema, std::string const &doc)
ResultKey const &	getPositiveKeys () const
	Return the standard flux keys registered by this algorithm. More...
ResultKey const &	getNegativeKeys () const
Public Member Functions inherited from lsst::meas::base::SimpleAlgorithm
virtual void	measureForced (afw::table::SourceRecord &measRecord, afw::image::Exposure< float > const &exposure, afw::table::SourceRecord const &refRecord, afw::image::Wcs const &refWcs) const
virtual void	measureNForced (afw::table::SourceCatalog const &measCat, afw::image::Exposure< float > const &exposure, afw::table::SourceCatalog const &refRecord, afw::image::Wcs const &refWcs) const
Public Member Functions inherited from lsst::meas::base::SingleFrameAlgorithm
virtual void	measureN (afw::table::SourceCatalog const &measCat, afw::image::Exposure< float > const &exposure) const
Public Member Functions inherited from lsst::meas::base::BaseAlgorithm
virtual	~BaseAlgorithm ()

Private Attributes
Control	_ctrl
afw::table::Key< float >	_chi2dofKey
meas::base::CentroidResultKey	_avgCentroid
meas::base::CentroidResultKey	_negCentroid
meas::base::CentroidResultKey	_posCentroid
afw::table::Key< afw::table::Flag >	_flagMaxPixelsKey

Additional Inherited Members
Protected Member Functions inherited from lsst::ip::diffim::DipoleFluxAlgorithm
	DipoleFluxAlgorithm (Control const &ctrl, std::string const &name, afw::table::Schema &schema, std::string const &doc, ResultKey const &positiveKeys, ResultKey const &negativeKeys)
	Initialize with a manually-constructed result key. More...
Protected Attributes inherited from lsst::ip::diffim::DipoleFluxAlgorithm
Control	_ctrl
meas::base::FluxResultKey	_fluxResultKey
meas::base::FlagHandler	_flagHandler
meas::base::SafeCentroidExtractor	_centroidExtractor
ResultKey	_positiveKeys
ResultKey	_negativeKeys

Detailed Description

Implementation of Psf dipole flux

Definition at line 314 of file DipoleAlgorithms.h.

Member Typedef Documentation

typedef PsfDipoleFluxControl lsst::ip::diffim::PsfDipoleFlux::Control

Definition at line 317 of file DipoleAlgorithms.h.

Constructor & Destructor Documentation

lsst::ip::diffim::PsfDipoleFlux::PsfDipoleFlux ( PsfDipoleFluxControl const &  ctrl, std::string const &  name, afw::table::Schema &  schema  )

inline

Definition at line 319 of file DipoleAlgorithms.h.

                                                                                                      :
        DipoleFluxAlgorithm(ctrl, name, schema, "jointly fitted psf flux counts"),
        _ctrl(ctrl),
        _chi2dofKey(schema.addField<float>(name+"_chi2dof",
                                           "chi2 per degree of freedom of fit")),
        _flagMaxPixelsKey(schema.addField<afw::table::Flag>(name+"_flags_maxpix",
                                                            "set if too large a footprint was sent to the algorithm"))
    {
        meas::base::CentroidResultKey::addFields(schema, name+"_pos_centroid", "psf fitted center of positive lobe", meas::base::SIGMA_ONLY);
        meas::base::CentroidResultKey::addFields(schema, name+"_neg_centroid", "psf fitted center of negative lobe", meas::base::SIGMA_ONLY);
        meas::base::CentroidResultKey::addFields(schema, name+"_centroid", "average of negative and positive lobe positions", meas::base::SIGMA_ONLY);
        _posCentroid = meas::base::CentroidResultKey(schema[name+"_pos_centroid"]);
        _negCentroid = meas::base::CentroidResultKey(schema[name+"_neg_centroid"]);
        _avgCentroid = meas::base::CentroidResultKey(schema[name+"_centroid"]);
    }

Member Function Documentation

std::pair< double, int > lsst::ip::diffim::PsfDipoleFlux::chi2	(	afw::table::SourceRecord &	source,
		afw::image::Exposure< float > const &	exposure,
		double	negCenterX,
		double	negCenterY,
		double	negFlux,
		double	posCenterX,
		double	poCenterY,
		double	posFlux
	)		const

Definition at line 295 of file DipoleAlgorithms.cc.

        {

    afw::geom::Point2D negCenter(negCenterX, negCenterY);
    afw::geom::Point2D posCenter(posCenterX, posCenterY);

    CONST_PTR(afw::detection::Footprint) footprint = source.getFootprint();

    /*
     * Fit for the superposition of Psfs at the two centroids.
     */
    CONST_PTR(afwDet::Psf) psf = exposure.getPsf();
    PTR(afwImage::Image<afwMath::Kernel::Pixel>) negPsf = psf->computeImage(negCenter);
    PTR(afwImage::Image<afwMath::Kernel::Pixel>) posPsf = psf->computeImage(posCenter);

    afwImage::Image<double> negModel(footprint->getBBox());
    afwImage::Image<double> posModel(footprint->getBBox());
    afwImage::Image<float> data(*(exposure.getMaskedImage().getImage()),footprint->getBBox());
    afwImage::Image<afwImage::VariancePixel> var(*(exposure.getMaskedImage().getVariance()),
                                                 footprint->getBBox());

    afwGeom::Box2I negPsfBBox = negPsf->getBBox();
    afwGeom::Box2I posPsfBBox = posPsf->getBBox();
    afwGeom::Box2I negModelBBox = negModel.getBBox();
    afwGeom::Box2I posModelBBox = posModel.getBBox();

    // Portion of the negative Psf that overlaps the model
    int negXmin = std::max(negPsfBBox.getMinX(), negModelBBox.getMinX());
    int negYmin = std::max(negPsfBBox.getMinY(), negModelBBox.getMinY());
    int negXmax = std::min(negPsfBBox.getMaxX(), negModelBBox.getMaxX());
    int negYmax = std::min(negPsfBBox.getMaxY(), negModelBBox.getMaxY());
    afwGeom::Box2I negBBox = afwGeom::Box2I(afwGeom::Point2I(negXmin, negYmin),
                                            afwGeom::Point2I(negXmax, negYmax));
    afwImage::Image<afwMath::Kernel::Pixel> negSubim(*negPsf, negBBox);
    afwImage::Image<double> negModelSubim(negModel, negBBox);
    negModelSubim += negSubim;

    // Portion of the positive Psf that overlaps the model
    int posXmin = std::max(posPsfBBox.getMinX(), posModelBBox.getMinX());
    int posYmin = std::max(posPsfBBox.getMinY(), posModelBBox.getMinY());
    int posXmax = std::min(posPsfBBox.getMaxX(), posModelBBox.getMaxX());
    int posYmax = std::min(posPsfBBox.getMaxY(), posModelBBox.getMaxY());
    afwGeom::Box2I posBBox = afwGeom::Box2I(afwGeom::Point2I(posXmin, posYmin),
                                            afwGeom::Point2I(posXmax, posYmax));
    afwImage::Image<afwMath::Kernel::Pixel> posSubim(*posPsf, posBBox);
    afwImage::Image<double> posModelSubim(posModel, posBBox);
    posModelSubim += posSubim;

    negModel  *= negFlux;   // scale negative model to image
    posModel  *= posFlux;   // scale positive model to image
    afwImage::Image<double> residuals(negModel, true); // full model contains negative lobe...
    residuals += posModel;  // plus positive lobe...
    residuals -= data;      // minus the data...
    residuals *= residuals; // squared...
    residuals /= var;       // divided by the variance : [(model-data)/sigma]**2
    afwMath::Statistics stats = afwMath::makeStatistics(residuals, afwMath::SUM | afwMath::NPOINT);
    double chi2 = stats.getValue(afwMath::SUM);
    int nPix = stats.getValue(afwMath::NPOINT);
    return std::pair<double,int>(chi2, nPix);
}

void lsst::ip::diffim::PsfDipoleFlux::fail ( afw::table::SourceRecord &  measRecord, 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 463 of file DipoleAlgorithms.cc.

                                                                                                  {
    _flagHandler.handleFailure(measRecord, error);
}

void lsst::ip::diffim::PsfDipoleFlux::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 360 of file DipoleAlgorithms.cc.

        {
    source.set(_flagMaxPixelsKey, true);

    typedef afw::image::Exposure<float>::MaskedImageT MaskedImageT;

    CONST_PTR(afw::detection::Footprint) footprint = source.getFootprint();
    if (!footprint) {
        throw LSST_EXCEPT(pex::exceptions::RuntimeError,
                          (boost::format("No footprint for source %d") % source.getId()).str());
    }

    if (footprint->getArea() > _ctrl.maxPixels) {
        // Too big
        return;
    }
    source.set(_flagMaxPixelsKey, false);

    afw::detection::PeakCatalog peakCatalog = afw::detection::PeakCatalog(footprint->getPeaks());

    if (peakCatalog.size() == 0) {
        throw LSST_EXCEPT(pex::exceptions::RuntimeError,
                          (boost::format("No peak for source %d") % source.getId()).str());
    }
    else if (peakCatalog.size() == 1) {
        // No deblending to do
        return;
    }

    // For N>=2, just measure the brightest-positive and brightest-negative
    // peaks.  peakCatalog is automatically ordered by peak flux, with the most
    // positive one (brightest) being first
    afw::detection::PeakRecord const& positivePeak = peakCatalog.front();
    afw::detection::PeakRecord const& negativePeak = peakCatalog.back();

    // Set up fit parameters and param names
    ROOT::Minuit2::MnUserParameters fitPar;

    fitPar.Add((boost::format("P%d")%NEGCENTXPAR).str(), negativePeak.getFx(), _ctrl.stepSizeCoord);
    fitPar.Add((boost::format("P%d")%NEGCENTYPAR).str(), negativePeak.getFy(), _ctrl.stepSizeCoord);
    fitPar.Add((boost::format("P%d")%NEGFLUXPAR).str(), negativePeak.getPeakValue(), _ctrl.stepSizeFlux);
    fitPar.Add((boost::format("P%d")%POSCENTXPAR).str(), positivePeak.getFx(), _ctrl.stepSizeCoord);
    fitPar.Add((boost::format("P%d")%POSCENTYPAR).str(), positivePeak.getFy(), _ctrl.stepSizeCoord);
    fitPar.Add((boost::format("P%d")%POSFLUXPAR).str(), positivePeak.getPeakValue(), _ctrl.stepSizeFlux);

    // Create the minuit object that knows how to minimise our functor
    //
    MinimizeDipoleChi2 minimizerFunc(*this, source, exposure);
    minimizerFunc.setErrorDef(_ctrl.errorDef);

    //
    // tell minuit about it
    //
    ROOT::Minuit2::MnMigrad migrad(minimizerFunc, fitPar);

    //
    // And let it loose
    //
    ROOT::Minuit2::FunctionMinimum min = migrad(_ctrl.maxFnCalls);

    float minChi2 = min.Fval();
    bool const isValid = min.IsValid() && std::isfinite(minChi2);

    if (true || isValid) {              // calculate coeffs even in minuit is unhappy

        /* I need to call chi2 one more time to grab nPix to calculate chi2/dof.
           Turns out that the Minuit operator method has to be const, and the
           measurement _apply method has to be const, so I can't store nPix as a
           private member variable anywhere.  Consted into a corner.
        */
        std::pair<double,int> fit = chi2(source, exposure,
                                         min.UserState().Value(NEGCENTXPAR),
                                         min.UserState().Value(NEGCENTYPAR),
                                         min.UserState().Value(NEGFLUXPAR),
                                         min.UserState().Value(POSCENTXPAR),
                                         min.UserState().Value(POSCENTYPAR),
                                         min.UserState().Value(POSFLUXPAR));
        double evalChi2 = fit.first;
        int nPix = fit.second;

        PTR(afw::geom::Point2D) minNegCentroid(new afw::geom::Point2D(min.UserState().Value(NEGCENTXPAR),
                                                                      min.UserState().Value(NEGCENTYPAR)));
        source.set(getNegativeKeys().getFlux(), min.UserState().Value(NEGFLUXPAR));
        source.set(getNegativeKeys().getFluxSigma(), min.UserState().Error(NEGFLUXPAR));

        PTR(afw::geom::Point2D) minPosCentroid(new afw::geom::Point2D(min.UserState().Value(POSCENTXPAR),
                                                                      min.UserState().Value(POSCENTYPAR)));
        source.set(getPositiveKeys().getFlux(), min.UserState().Value(POSFLUXPAR));
        source.set(getPositiveKeys().getFluxSigma(), min.UserState().Error(POSFLUXPAR));

        source.set(_chi2dofKey, evalChi2 / (nPix - minimizerFunc.getNpar()));
        source.set(_negCentroid.getX(), minNegCentroid->getX());
        source.set(_negCentroid.getY(), minNegCentroid->getY());
        source.set(_posCentroid.getX(), minPosCentroid->getX());
        source.set(_posCentroid.getY(), minPosCentroid->getY());
        source.set(_avgCentroid.getX(), 0.5*(minNegCentroid->getX() + minPosCentroid->getX()));
        source.set(_avgCentroid.getY(), 0.5*(minNegCentroid->getY() + minPosCentroid->getY()));

    }
}

Member Data Documentation

meas::base::CentroidResultKey lsst::ip::diffim::PsfDipoleFlux::_avgCentroid

private

Definition at line 354 of file DipoleAlgorithms.h.

afw::table::Key<float> lsst::ip::diffim::PsfDipoleFlux::_chi2dofKey

private

Definition at line 353 of file DipoleAlgorithms.h.

Control lsst::ip::diffim::PsfDipoleFlux::_ctrl

private

Definition at line 352 of file DipoleAlgorithms.h.

afw::table::Key< afw::table::Flag > lsst::ip::diffim::PsfDipoleFlux::_flagMaxPixelsKey

private

Definition at line 358 of file DipoleAlgorithms.h.

meas::base::CentroidResultKey lsst::ip::diffim::PsfDipoleFlux::_negCentroid

private

Definition at line 355 of file DipoleAlgorithms.h.

meas::base::CentroidResultKey lsst::ip::diffim::PsfDipoleFlux::_posCentroid

private

Definition at line 356 of file DipoleAlgorithms.h.

---

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

• /home/lsstsw/stack/Linux64/ip_diffim/11.0+1/include/lsst/ip/diffim/DipoleAlgorithms.h
• /home/lsstsw/stack/Linux64/ip_diffim/11.0+1/src/DipoleAlgorithms.cc