GeneralPsfFitter.h

Go to the documentation of this file.

// -*- lsst-c++ -*-
/*
 * LSST Data Management System
 * Copyright 2008-2013 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_MODELFIT_GeneralPsfFitter_h_INCLUDED
#define LSST_MEAS_MODELFIT_GeneralPsfFitter_h_INCLUDED
 
#include <memory>
 
#include "lsst/pex/config.h"
#include "lsst/shapelet/FunctorKeys.h"
#include "lsst/meas/modelfit/Model.h"
#include "lsst/meas/modelfit/Prior.h"
#include "lsst/meas/modelfit/Likelihood.h"
#include "lsst/geom.h"
#include "lsst/afw/geom.h"
#include "lsst/afw/table/Source.h"
#include "lsst/meas/base/exceptions.h"
#include "lsst/meas/base/FlagHandler.h"
#include "lsst/meas/modelfit/optimizer.h"
 
namespace lsst { namespace meas { namespace modelfit {

class GeneralPsfFitterComponentControl {
public:
 
    GeneralPsfFitterComponentControl(int order_=0, double radiusFactor_=1.0) :
        order(order_), positionPriorSigma(0.1), ellipticityPriorSigma(0.3),
        radiusFactor(radiusFactor_), radiusPriorSigma(0.5)
    {}
 
    LSST_CONTROL_FIELD(
        order, int,
        "shapelet order for this component; negative to disable this component completely"
    );
    LSST_CONTROL_FIELD(
        positionPriorSigma, double,
        "sigma (in pixels) in an isotropic 2-d Gaussian prior on the center of this shapelet component, "
        "relative to the center of the PSF image"
    );
    LSST_CONTROL_FIELD(
        ellipticityPriorSigma, double,
        "sigma in an isotropic 2-d Gaussian prior on the conformal-shear ellipticity eta"
    );
    LSST_CONTROL_FIELD(
        radiusFactor, double,
        "Sets the fiducial radius of this component relative to the 'primary radius' of the PSF: either "
        "the second-moments radius of the PSF image (in an initial fit), or the radius of the primary "
        "component in a previous fit.  Ignored if the previous fit included this component (as then we "
        "can just use that radius)."
    );
    LSST_CONTROL_FIELD(
        radiusPriorSigma, double,
        "sigma in a Gaussian prior on ln(radius/fiducialRadius)"
    );
 
};

class GeneralPsfFitterControl {
public:
 
    GeneralPsfFitterControl() :
        inner(-1, 0.5), primary(0, 1.0), wings(0, 2.0), outer(-1, 4.0), defaultNoiseSigma(0.001)
    {}
 
    LSST_NESTED_CONTROL_FIELD(
        inner, lsst.meas.modelfit.psf.psfContinued, GeneralPsfFitterComponentControl,
        "Innermost shapelet expansion, used to fit PSFs with very sharp cores"
    );
 
    LSST_NESTED_CONTROL_FIELD(
        primary, lsst.meas.modelfit.psf.psfContinued, GeneralPsfFitterComponentControl,
        "Primary shapelet expansion, typically used to fit the bulk of the PSF "
    );
 
    LSST_NESTED_CONTROL_FIELD(
        wings, lsst.meas.modelfit.psf.psfContinued, GeneralPsfFitterComponentControl,
        "Wing shapelet expansion (between primary and outer), typically used to fit the wings of the PSF"
    );
 
    LSST_NESTED_CONTROL_FIELD(
        outer, lsst.meas.modelfit.psf.psfContinued, GeneralPsfFitterComponentControl,
        "Outermost shapelet expansion, used to fit PSFs with very broad wings"
    );
 
    LSST_NESTED_CONTROL_FIELD(
        optimizer, lsst.meas.modelfit.optimizer, OptimizerControl,
        "Configuration of the optimizer used to do the fitting"
    );
 
    LSST_CONTROL_FIELD(
        defaultNoiseSigma, double, "Default value for the noiseSigma parameter in GeneralPsfFitter.apply()"
    );
 
};

class GeneralPsfFitter {
public:
    explicit GeneralPsfFitter(GeneralPsfFitterControl const & ctrl);

    shapelet::MultiShapeletFunctionKey addFields(
        afw::table::Schema & schema,
        std::string const & prefix
    ) const;

    std::shared_ptr<Model> getModel() const { return _model; }

    std::shared_ptr<Prior> getPrior() const { return _prior; }

    shapelet::MultiShapeletFunction adapt(
        shapelet::MultiShapeletFunction const & previousFit,
        std::shared_ptr<Model> previousModel
    ) const;


    shapelet::MultiShapeletFunction apply(
        afw::image::Image<Pixel> const & image,
        afw::geom::ellipses::Quadrupole const & moments,
        Scalar noiseSigma=-1,
        int * pState = nullptr
    ) const;
    shapelet::MultiShapeletFunction apply (
        afw::image::Image<double> const & image,
        afw::geom::ellipses::Quadrupole const & moments,
        Scalar noiseSigma=-1,
        int * pState = nullptr
    ) const {
        return apply(afw::image::Image<float>(image, true), moments, noiseSigma, pState);
    }



    shapelet::MultiShapeletFunction apply(
        afw::image::Image<Pixel> const & image,
        shapelet::MultiShapeletFunction const & initial,
        Scalar noiseSigma=-1,
        int * pState = nullptr
    ) const;
    shapelet::MultiShapeletFunction apply(
        afw::image::Image<double> const & image,
        shapelet::MultiShapeletFunction const & initial,
        Scalar noiseSigma=-1,
        int * pState = nullptr
    ) const {
        return apply(afw::image::Image<float>(image, true), initial, noiseSigma);
    }

 
private:
    GeneralPsfFitterControl _ctrl;
    std::shared_ptr<Model> _model;
    std::shared_ptr<Prior> _prior;
};
 
class GeneralPsfFitterAlgorithm : public GeneralPsfFitter {
public:
 
    // Structures and routines to manage flaghandler
    static base::FlagDefinitionList const & getFlagDefinitions();
    static base::FlagDefinition const FAILURE;
    static base::FlagDefinition const MAX_INNER_ITERATIONS;
    static base::FlagDefinition const MAX_OUTER_ITERATIONS;
    static base::FlagDefinition const EXCEPTION;
    static base::FlagDefinition const CONTAINS_NAN;
 
    typedef GeneralPsfFitterControl Control;
 
    GeneralPsfFitterAlgorithm(GeneralPsfFitterControl const & ctrl,
        afw::table::Schema & schema,
        std::string const & prefix
    );
 
    shapelet::MultiShapeletFunctionKey getKey() {
        return _key;
    }
 
    void measure(
        afw::table::SourceRecord & measRecord,
        afw::image::Image<double> const & image,
        shapelet::MultiShapeletFunction const & initial
    ) const;
 
    void measure(
        afw::table::SourceRecord & measRecord,
        afw::image::Image<double> const & image,
        afw::geom::ellipses::Quadrupole const & moments
    ) const;
 
    void fail(
        afw::table::SourceRecord & measRecord,
        lsst::meas::base::MeasurementError * error=nullptr
    ) const;
 
private:
    shapelet::MultiShapeletFunctionKey _key;
    lsst::meas::base::FlagHandler _flagHandler;
};

class MultiShapeletPsfLikelihood : public Likelihood {
public:
 
    MultiShapeletPsfLikelihood(
        ndarray::Array<Pixel const,2,1> const & image,
        geom::Point2I const & xy0,
        std::shared_ptr<Model> model,
        Scalar sigma,
        ndarray::Array<Scalar const,1,1> const & fixed
    );
 
    void computeModelMatrix(
        ndarray::Array<Pixel,2,-1> const & modelMatrix,
        ndarray::Array<Scalar const,1,1> const & nonlinear,
        bool doApplyWeights=true
    ) const override;
 
    virtual ~MultiShapeletPsfLikelihood();
 
private:
    class Impl;
    std::unique_ptr<Impl> _impl;
};
 
}}} // namespace lsst::meas::modelfit
 
#endif // !LSST_MEAS_MODELFIT_GeneralPsfFitter_h_INCLUDED