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;

    PTR(Model) getModel() const { return _model; }

    PTR(Prior) getPrior() const { return _prior; }

    shapelet::MultiShapeletFunction adapt(
        shapelet::MultiShapeletFunction const & previousFit,
        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;
    PTR(Model) _model;
    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,
        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