photometryKron.h

Go to the documentation of this file.

// -*- lsst-c++ -*-
/*
 * LSST Data Management System
 * Copyright 2008-2015 LSST Corporation.
 *
 * This product includes software developed by the
 * LSST Project (http://www.lsst.org/).
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the LSST License Statement and
 * the GNU General Public License along with this program.  If not,
 * see <http://www.lsstcorp.org/LegalNotices/>.
 */
#ifndef LSST_MEAS_EXTENSIONS_PHOTOMETRY_KRON_H
#define LSST_MEAS_EXTENSIONS_PHOTOMETRY_KRON_H
 
#include <memory>
#include <cmath>
 
#include "lsst/pex/config.h"
#include "lsst/geom.h"
#include "lsst/afw/image/Exposure.h"
#include "lsst/afw/geom/SkyWcs.h"
#include "lsst/meas/base/Algorithm.h"
#include "lsst/meas/base/FluxUtilities.h"
#include "lsst/meas/base/CentroidUtilities.h"
#include "lsst/meas/base/FlagHandler.h"
#include "lsst/meas/base/InputUtilities.h"
 
namespace lsst { namespace meas { namespace extensions { namespace photometryKron {
 
class KronAperture;
 
class KronFluxControl {
public:
 
    LSST_CONTROL_FIELD(fixed, bool,
                       "if true, use existing shape and centroid measurements instead of fitting");
    LSST_CONTROL_FIELD(nSigmaForRadius, double,
                       "Multiplier of rms size for aperture used to initially estimate the Kron radius");
    LSST_CONTROL_FIELD(nIterForRadius, int, "Number of times to iterate when setting the Kron radius");
    LSST_CONTROL_FIELD(nRadiusForFlux, double, "Number of Kron radii for Kron flux");
    LSST_CONTROL_FIELD(maxSincRadius, double,
                       "Largest aperture for which to use the slow, accurate, sinc aperture code");
    LSST_CONTROL_FIELD(minimumRadius, double,
                       "Minimum Kron radius (if == 0.0 use PSF's Kron radius) if enforceMinimumRadius. "
                       "Also functions as fallback aperture radius if set.");
    LSST_CONTROL_FIELD(enforceMinimumRadius, bool, "If true check that the Kron radius exceeds some minimum");
    LSST_CONTROL_FIELD(useFootprintRadius, bool,
                       "Use the Footprint size as part of initial estimate of Kron radius");
    LSST_CONTROL_FIELD(smoothingSigma, double,
                       "Smooth image with N(0, smoothingSigma^2) Gaussian while estimating R_K");
    LSST_CONTROL_FIELD(refRadiusName, std::string,
                       "Name of field specifying reference Kron radius for forced measurement");
 
    KronFluxControl() :
        fixed(false),
        nSigmaForRadius(6.0),
        nIterForRadius(1),
        nRadiusForFlux(2.5),
        maxSincRadius(10.0),
        minimumRadius(0.0),
        enforceMinimumRadius(true),
        useFootprintRadius(false),
        smoothingSigma(-1.0),
        refRadiusName("ext_photometryKron_KronFlux_radius")
    {}
};
 
class KronFluxAlgorithm : public base::SimpleAlgorithm {
public:
 
    // Structures and routines to manage flaghandler
    static meas::base::FlagDefinitionList const & getFlagDefinitions();
    static meas::base::FlagDefinition const FAILURE;
    static meas::base::FlagDefinition const EDGE;
    static meas::base::FlagDefinition const BAD_SHAPE_NO_PSF;
    static meas::base::FlagDefinition const NO_MINIMUM_RADIUS;
    static meas::base::FlagDefinition const NO_FALLBACK_RADIUS;
    static meas::base::FlagDefinition const BAD_RADIUS;
    static meas::base::FlagDefinition const USED_MINIMUM_RADIUS;
    static meas::base::FlagDefinition const USED_PSF_RADIUS;
    static meas::base::FlagDefinition const SMALL_RADIUS;
    static meas::base::FlagDefinition const BAD_SHAPE;
 
    typedef KronFluxControl Control;
 
    KronFluxAlgorithm(
        Control const & ctrl,
        std::string const & name,
        afw::table::Schema & schema,
        daf::base::PropertySet & metadata
    );
 
    virtual void measure(
        afw::table::SourceRecord & measRecord,
        afw::image::Exposure<float> const & exposure
    ) const;
 
    virtual void measureForced(
        afw::table::SourceRecord & measRecord,
        afw::image::Exposure<float> const & exposure,
        afw::table::SourceRecord const & refRecord,
        afw::geom::SkyWcs const & refWcs
    ) const;
 
    virtual void fail(
        afw::table::SourceRecord & measRecord,
        meas::base::MeasurementError * error=NULL
    ) const;
 
private:
 
    void _applyAperture(
        afw::table::SourceRecord & source,
        afw::image::Exposure<float> const& exposure,
        KronAperture const& aperture
        ) const;
 
    void _applyForced(
        afw::table::SourceRecord & source,
        afw::image::Exposure<float> const & exposure,
        geom::Point2D const & center,
        afw::table::SourceRecord const & reference,
        geom::AffineTransform const & refToMeas
    ) const;
 
    std::shared_ptr<KronAperture> _fallbackRadius(afw::table::SourceRecord& source, double const R_K_psf,
                                      pex::exceptions::Exception& exc) const;
 
    std::string _name;
    Control _ctrl;
    meas::base::FluxResultKey _fluxResultKey;
    afw::table::Key<float> _radiusKey;
    afw::table::Key<float> _radiusForRadiusKey;
    afw::table::Key<float> _psfRadiusKey;
    meas::base::FlagHandler _flagHandler;
    meas::base::SafeCentroidExtractor _centroidExtractor;
};
 
class KronAperture {
public:
    KronAperture(geom::Point2D const& center, afw::geom::ellipses::BaseCore const& core,
                 float radiusForRadius=std::nanf("")) :
        _center(center),
        _axes(core),
        _radiusForRadius(radiusForRadius)
        {}
 
    explicit KronAperture(afw::table::SourceRecord const& source, float radiusForRadius=std::nanf("")) :
        _center(geom::Point2D(source.getX(), source.getY())),
        _axes(source.getShape()),
        _radiusForRadius(radiusForRadius)
        {}
 
    KronAperture(afw::table::SourceRecord const& reference, geom::AffineTransform const& refToMeas,
                 double radius, float radiusForRadius=std::nanf("")) :
        _center(refToMeas(reference.getCentroid())),
        _axes(getKronAxes(reference.getShape(), refToMeas.getLinear(), radius)),
        _radiusForRadius(radiusForRadius)
        {}
 
    double getX() const { return _center.getX(); }
    double getY() const { return _center.getY(); }
    float getRadiusForRadius() const { return _radiusForRadius; }
 
    geom::Point2D const& getCenter() const { return _center; }
 
    afw::geom::ellipses::Axes & getAxes() { return _axes; }
 
    afw::geom::ellipses::Axes const& getAxes() const { return _axes; }
 
 
    template<typename ImageT>
    static std::shared_ptr<KronAperture> determineRadius(
        ImageT const& image,  
        afw::geom::ellipses::Axes axes,  
        geom::Point2D const& center,   
        KronFluxControl const& ctrl  
        );
 
    template<typename ImageT>
    std::pair<double, double> measureFlux(
        ImageT const& image,  
        double const nRadiusForFlux,  
        double const maxSincRadius  
        ) const;
 
    std::shared_ptr<KronAperture> transform(geom::AffineTransform const& trans) const {
        geom::Point2D const center = trans(getCenter());
        afw::geom::ellipses::Axes const axes(*getAxes().transform(trans.getLinear()).copy());
        return std::make_shared<KronAperture>(center, axes);
    }
 
    static
    afw::geom::ellipses::Axes getKronAxes(
        afw::geom::ellipses::Axes const& shape,
        geom::LinearTransform const& transformation,
        double const radius
        );
 
private:
    geom::Point2D const _center;     // Center of aperture
    afw::geom::ellipses::Axes _axes;      // Ellipse defining aperture shape
    float _radiusForRadius;               // Radius used to estimate the Kron radius
};
 
}}}} // namespace lsst::meas::extensions::photometryKron
 
#endif // !LSST_MEAS_EXTENSIONS_PHOTOMETRY_KRON_H