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

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

Public Member Functions
	PsfDipoleFlux (PsfDipoleFluxControl const &ctrl, afw::table::Schema &schema)
template<typename PixelT >
std::pair< double, int >	chi2 (afw::table::SourceRecord &source, afw::image::Exposure< PixelT > const &exposure, double negCenterX, double negCenterY, double negFlux, double posCenterX, double poCenterY, double posFlux) const
Public Member Functions inherited from lsst::ip::diffim::DipoleFluxAlgorithm
DipoleFluxControl const &	getControl () const
	Return a clone of the control object used to construct the algorithm. More...
KeyTuple const &	getPositiveKeys () const
	Return the standard flux keys registered by this algorithm. More...
KeyTuple const &	getNegativeKeys () const
Public Member Functions inherited from lsst::meas::algorithms::Algorithm
	Algorithm (AlgorithmControl const &ctrl)
virtual	~Algorithm ()
AlgorithmControl const &	getControl () const
	Return a clone of the control object used to construct the algorithm. More...
template<typename PixelT >
void	apply (afw::table::SourceRecord &source, afw::image::Exposure< PixelT > const &exposure, afw::geom::Point2D const &center) const
	Run the algorithm, filling appropriate fields in the given source. More...

Private Member Functions
template<typename PixelT >
void	_apply (afw::table::SourceRecord &source, afw::image::Exposure< PixelT > const &exposure, afw::geom::Point2D const &center) const
	LSST_MEAS_ALGORITHM_PRIVATE_INTERFACE (PsfDipoleFlux)

Private Attributes
afw::table::Key< float >	_chi2dofKey
afw::table::KeyTuple < afw::table::Centroid >	_avgCentroid
afw::table::KeyTuple < afw::table::Centroid >	_negCentroid
afw::table::KeyTuple < afw::table::Centroid >	_posCentroid
afw::table::Key< afw::table::Flag >	_flagMaxPixelsKey

Additional Inherited Members
Public Types inherited from lsst::ip::diffim::DipoleFluxAlgorithm
typedef afw::table::KeyTuple < afw::table::Flux >	KeyTuple
	Tuple type that holds the keys that define a standard flux algorithm. More...
Protected Member Functions inherited from lsst::ip::diffim::DipoleFluxAlgorithm
	DipoleFluxAlgorithm (DipoleFluxControl const &ctrl, KeyTuple const &positiveKeys, KeyTuple const &negativeKeys)
	Initialize with a manually-constructed key tuple. More...
	DipoleFluxAlgorithm (DipoleFluxControl const &ctrl, afw::table::Schema &schema, char const *doc)
	Initialize using afw::table::addFlux field to fill out repetitive descriptions. More...
Protected Member Functions inherited from lsst::meas::algorithms::Algorithm
template<typename PixelT >
void	_apply (afw::table::SourceRecord &source, afw::image::Exposure< PixelT > const &exposure, afw::geom::Point2D const &center) const
	Simulated virtual function that all algorithms must implement. More...

Detailed Description

Implementation of Psf dipole flux

Definition at line 324 of file dipoleAlgorithm.cc.

Constructor & Destructor Documentation

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

inline

Definition at line 327 of file dipoleAlgorithm.cc.

                                                                              :
        DipoleFluxAlgorithm(ctrl, schema, "jointly fitted psf flux counts"),
        _chi2dofKey(schema.addField<float>(ctrl.name+".chi2dof", 
                                           "chi2 per degree of freedom of fit")),
        _avgCentroid(
            addCentroidFields(schema, ctrl.name+".centroid", 
                              "average of the postive and negative lobe positions")),
        _negCentroid(
            addCentroidFields(schema, ctrl.name+".neg.centroid", 
                              "psf fitted center of negative lobe")),
        _posCentroid(
            addCentroidFields(schema, ctrl.name+".pos.centroid", 
                              "psf fitted center of positive lobe")),
        _flagMaxPixelsKey(schema.addField<afw::table::Flag>(ctrl.name+".flags.maxpix",
                                                            "set if too large a footprint was sent to the algorithm"))
    {}

Member Function Documentation

template<typename PixelT >

void lsst::ip::diffim::PsfDipoleFlux::_apply ( afw::table::SourceRecord &  source, afw::image::Exposure< PixelT > const &  exposure, afw::geom::Point2D const &  center  ) const

private

Definition at line 495 of file dipoleAlgorithm.cc.

        {

    source.set(getPositiveKeys().flag, true); // say we've failed so that's the result if we throw
    source.set(getNegativeKeys().flag, true); // say we've failed so that's the result if we throw
    source.set(_flagMaxPixelsKey, true);

    typedef typename afw::image::Exposure<PixelT>::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());
    }

    PsfDipoleFluxControl const & ctrl = static_cast<PsfDipoleFluxControl const &>(getControl());
    if (footprint->getArea() > ctrl.maxPixels) {
        // Too big
        return;
    }
    source.set(_flagMaxPixelsKey, false);

    afw::detection::Footprint::PeakList peakList = 
        afw::detection::Footprint::PeakList(footprint->getPeaks());

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

    // For N>=2, just measure the brightest-positive and brightest-negative
    // peaks.  peakList is automatically ordered by peak flux, with the most
    // positive one (brightest) being first
    PTR(afw::detection::Peak) positivePeak = peakList.front();
    PTR(afw::detection::Peak) negativePeak = peakList.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<PixelT> 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().meas, min.UserState().Value(NEGFLUXPAR));
        source.set(getNegativeKeys().err, min.UserState().Error(NEGFLUXPAR));
        source.set(getNegativeKeys().flag, false);
        
        PTR(afw::geom::Point2D) minPosCentroid(new afw::geom::Point2D(min.UserState().Value(POSCENTXPAR), 
                                                                      min.UserState().Value(POSCENTYPAR))); 
        source.set(getPositiveKeys().meas, min.UserState().Value(POSFLUXPAR));
        source.set(getPositiveKeys().err, min.UserState().Error(POSFLUXPAR));
        source.set(getPositiveKeys().flag, false);

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

    }

}

template<typename PixelT >

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

Definition at line 428 of file dipoleAlgorithm.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<PixelT> 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);
}    

lsst::ip::diffim::PsfDipoleFlux::LSST_MEAS_ALGORITHM_PRIVATE_INTERFACE ( PsfDipoleFlux  )

private

Member Data Documentation

afw::table::KeyTuple<afw::table::Centroid> lsst::ip::diffim::PsfDipoleFlux::_avgCentroid

private

Definition at line 362 of file dipoleAlgorithm.cc.

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

private

Definition at line 361 of file dipoleAlgorithm.cc.

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

private

Definition at line 366 of file dipoleAlgorithm.cc.

afw::table::KeyTuple<afw::table::Centroid> lsst::ip::diffim::PsfDipoleFlux::_negCentroid

private

Definition at line 363 of file dipoleAlgorithm.cc.

afw::table::KeyTuple<afw::table::Centroid> lsst::ip::diffim::PsfDipoleFlux::_posCentroid

private

Definition at line 364 of file dipoleAlgorithm.cc.

---

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

• /lsst/home/lsstsw/stack/Linux64/ip_diffim/master-gcec1da163f+63/src/dipoleAlgorithm.cc