ShapeletFunction.h

Go to the documentation of this file.

// -*- LSST-C++ -*-

/*
 * LSST Data Management System
 * Copyright 2008, 2009, 2010, 2011 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_AFW_MATH_SHAPELETS_SHAPELETFUNCTION_H
#define LSST_AFW_MATH_SHAPELETS_SHAPELETFUNCTION_H

#include "ndarray.h"
#include "lsst/shapelet/constants.h"
#include "lsst/shapelet/GaussHermiteEvaluator.h"
#include "lsst/shapelet/ConversionMatrix.h"
#include "lsst/geom.h"
#include "lsst/afw/geom/ellipses.h"
#include "lsst/afw/image/Image.h"

#include <list>

namespace lsst { namespace shapelet {

class ShapeletFunctionEvaluator;

class ShapeletFunction {
public:

    typedef std::shared_ptr<ShapeletFunction> Ptr;
    typedef std::shared_ptr<ShapeletFunction const> ConstPtr;

    typedef ShapeletFunctionEvaluator Evaluator;

    static double const FLUX_FACTOR;

    int getOrder() const { return _order; }

    afw::geom::ellipses::Ellipse const & getEllipse() const { return _ellipse; }

    afw::geom::ellipses::Ellipse & getEllipse() { return _ellipse; }

    void setEllipse(afw::geom::ellipses::Ellipse const & ellipse) { _ellipse = ellipse; }

    BasisTypeEnum getBasisType() const { return _basisType; }

    void changeBasisType(BasisTypeEnum basisType) {
        ConversionMatrix::convertCoefficientVector(_coefficients, _basisType, basisType, _order);
        _basisType = basisType;
    }

    void normalize(double value=1.0);

    ndarray::Array<double,1,1> const getCoefficients() { return _coefficients; }

    ndarray::Array<double const,1,1> const getCoefficients() const { return _coefficients; }

    ShapeletFunction convolve(ShapeletFunction const & other) const;

    Evaluator evaluate() const;

    void shiftInPlace(geom::Extent2D const & offset) {
        _ellipse.getCenter() += offset;
    }

    void transformInPlace(geom::AffineTransform const & transform) {
        _ellipse.transform(transform).inPlace();
    }

    ShapeletFunction(int order, BasisTypeEnum basisType);

    ShapeletFunction(
        int order, BasisTypeEnum basisType,
        ndarray::Array<double,1,1> const & coefficients
    );

    ShapeletFunction(int order, BasisTypeEnum basisType, double radius,
                     geom::Point2D const & center=geom::Point2D());

    ShapeletFunction(
        int order, BasisTypeEnum basisType, double radius, geom::Point2D const & center,
        ndarray::Array<double,1,1> const & coefficients
    );

    ShapeletFunction(int order, BasisTypeEnum basisType,
        afw::geom::ellipses::Ellipse const & ellipse);

    ShapeletFunction(
        int order, BasisTypeEnum basisType,
        afw::geom::ellipses::Ellipse const & ellipse,
        ndarray::Array<double const,1,1> const & coefficients
    );

    ShapeletFunction(ShapeletFunction const & other);

    ShapeletFunction & operator=(ShapeletFunction const & other);

    ShapeletFunction();

private:

    int _order;
    BasisTypeEnum _basisType;
    lsst::afw::geom::ellipses::Ellipse _ellipse;
    ndarray::Array<double,1,1> _coefficients;
};

class ShapeletFunctionEvaluator {
public:

    typedef std::shared_ptr<ShapeletFunctionEvaluator> Ptr;
    typedef std::shared_ptr<ShapeletFunctionEvaluator const> ConstPtr;

    double operator()(double x, double y) const {
        return this->operator()(geom::Point2D(x, y));
    }

    double operator()(geom::Point2D const & point) const {
        return _normalization * _h.sumEvaluation(_coefficients, _transform(point));
    }

    double operator()(geom::Extent2D const & point) const {
        return _normalization * _h.sumEvaluation(_coefficients, _transform(point));
    }

    ndarray::Array<double,1,1> operator()(
        ndarray::Array<double const,1> const & x,
        ndarray::Array<double const,1> const & y
    ) const;

    void addToImage(
        ndarray::Array<double,2,1> const & array,
        geom::Point2I const & xy0 = geom::Point2I()
    ) const;

    void addToImage(afw::image::Image<double> & image) const {
        addToImage(image.getArray(), image.getXY0());
    }

    double integrate() const {
        return _h.sumIntegration(_coefficients);
    }

    afw::geom::ellipses::Ellipse computeMoments() const;

    void update(ShapeletFunction const & function);

    explicit ShapeletFunctionEvaluator(ShapeletFunction const & function);

private:

    friend class MultiShapeletFunctionEvaluator;

    void _computeRawMoments(double & q0, Eigen::Vector2d & q1, Eigen::Matrix2d & q2) const;

    double _normalization;
    ndarray::Array<double const,1,1> _coefficients;
    geom::AffineTransform _transform;
    GaussHermiteEvaluator _h;
};

inline ShapeletFunctionEvaluator ShapeletFunction::evaluate() const {
    return ShapeletFunctionEvaluator(*this);
}

}} // namespace lsst::shapelet

#endif // !defined(LSST_AFW_MATH_SHAPELETS_SHAPELETFUNCTION_H)