Kernel.h

Go to the documentation of this file.

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

/*
 * LSST Data Management System
 * Copyright 2008-2016  AURA/LSST.
 *
 * 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_KERNEL_H
#define LSST_AFW_MATH_KERNEL_H
/*
 * Declare the Kernel class and subclasses.
 */
#include <memory>
#include <type_traits>
#include <utility>
#include <vector>

#include "boost/mpl/or.hpp"

#include "lsst/daf/base/Citizen.h"
#include "lsst/geom.h"
#include "lsst/afw/image/Image.h"
#include "lsst/afw/image/Utils.h"
#include "lsst/afw/math/Function.h"
#include "lsst/afw/math/traits.h"

#include "lsst/afw/table/io/Persistable.h"

namespace lsst {
namespace afw {

namespace math {

class Kernel : public lsst::daf::base::Citizen,
               public afw::table::io::PersistableFacade<Kernel>,
               public afw::table::io::Persistable {
public:
    typedef double Pixel;
    typedef std::shared_ptr<lsst::afw::math::Function2<double>> SpatialFunctionPtr;
    typedef lsst::afw::math::Function2<double> SpatialFunction;
    typedef lsst::afw::math::NullFunction2<double> NullSpatialFunction;

    // Traits values for this class of Kernel
    typedef generic_kernel_tag kernel_fill_factor;

    explicit Kernel();

    explicit Kernel(int width, int height, unsigned int nKernelParams,
                    SpatialFunction const &spatialFunction = NullSpatialFunction());

    explicit Kernel(int width, int height, const std::vector<SpatialFunctionPtr> spatialFunctionList);

    // prevent copying and assignment (to avoid problems from type slicing)
    Kernel(const Kernel &) = delete;
    Kernel(Kernel &&) = delete;
    Kernel &operator=(const Kernel &) = delete;
    Kernel &operator=(Kernel &&) = delete;

    ~Kernel() override = default;

    virtual std::shared_ptr<Kernel> clone() const = 0;

    virtual std::shared_ptr<Kernel> resized(int width, int height) const = 0;

    double computeImage(lsst::afw::image::Image<Pixel> &image, bool doNormalize, double x = 0.0,
                        double y = 0.0) const;

    lsst::geom::Extent2I const getDimensions() const { return lsst::geom::Extent2I(_width, _height); }

    void setDimensions(lsst::geom::Extent2I dims) {
        _width = dims.getX();
        _height = dims.getY();
    }
    inline void setWidth(int width) { _width = width; }
    inline void setHeight(int height) { _height = height; }

    inline int getWidth() const { return _width; }

    inline int getHeight() const { return _height; }

    inline lsst::geom::Point2I getCtr() const { return lsst::geom::Point2I(_ctrX, _ctrY); }

    inline int getCtrX() const { return _ctrX; }

    inline int getCtrY() const { return _ctrY; }

    inline lsst::geom::Box2I getBBox() const {
        return lsst::geom::Box2I(lsst::geom::Point2I(-_ctrX, -_ctrY), lsst::geom::Extent2I(_width, _height));
    }

    inline unsigned int getNKernelParameters() const { return _nKernelParams; }

    inline int getNSpatialParameters() const {
        return this->isSpatiallyVarying() ? _spatialFunctionList[0]->getNParameters() : 0;
    }

    SpatialFunctionPtr getSpatialFunction(unsigned int index) const;

    std::vector<SpatialFunctionPtr> getSpatialFunctionList() const;

    virtual double getKernelParameter(unsigned int i) const { return getKernelParameters()[i]; }

    virtual std::vector<double> getKernelParameters() const;

    lsst::geom::Box2I growBBox(lsst::geom::Box2I const &bbox) const;

    lsst::geom::Box2I shrinkBBox(lsst::geom::Box2I const &bbox) const;

    inline void setCtr(lsst::geom::Point2I ctr) {
        _ctrX = ctr.getX();
        _ctrY = ctr.getY();
        _setKernelXY();
    }

    inline void setCtrX(int ctrX) {
        _ctrX = ctrX;
        _setKernelXY();
    }

    inline void setCtrY(int ctrY) {
        _ctrY = ctrY;
        _setKernelXY();
    }

    inline std::vector<std::vector<double>> getSpatialParameters() const {
        std::vector<std::vector<double>> spatialParams;
        std::vector<SpatialFunctionPtr>::const_iterator spFuncIter = _spatialFunctionList.begin();
        for (; spFuncIter != _spatialFunctionList.end(); ++spFuncIter) {
            spatialParams.push_back((*spFuncIter)->getParameters());
        }
        return spatialParams;
    }

    inline bool isSpatiallyVarying() const { return _spatialFunctionList.size() != 0; }

    inline void setKernelParameters(std::vector<double> const &params) {
        if (this->isSpatiallyVarying()) {
            throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError, "Kernel is spatially varying");
        }
        const unsigned int nParams = this->getNKernelParameters();
        if (nParams != params.size()) {
            throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterError,
                              (boost::format("Number of parameters is wrong, saw %d expected %d") % nParams %
                               params.size())
                                      .str());
        }
        for (unsigned int ii = 0; ii < nParams; ++ii) {
            this->setKernelParameter(ii, params[ii]);
        }
    }

    inline void setKernelParameters(std::pair<double, double> const &params) {
        this->setKernelParameter(0, params.first);
        this->setKernelParameter(1, params.second);
    }

    void setSpatialParameters(const std::vector<std::vector<double>> params);

    void computeKernelParametersFromSpatialModel(std::vector<double> &kernelParams, double x, double y) const;

    virtual std::string toString(std::string const &prefix = "") const;

    virtual void computeCache(int const  
    ) {}

    virtual int getCacheSize() const { return 0; };

#if 0  // fails to compile with icc; is it actually used?
        virtual void toFile(std::string fileName) const;
#endif

    struct PersistenceHelper;

protected:
    std::string getPythonModule() const override;

    virtual void setKernelParameter(unsigned int ind, double value) const;

    void setKernelParametersFromSpatialModel(double x, double y) const;

    virtual double doComputeImage(lsst::afw::image::Image<Pixel> &image, bool doNormalize) const = 0;

    std::vector<SpatialFunctionPtr> _spatialFunctionList;

private:
    int _width;
    int _height;
    int _ctrX;
    int _ctrY;
    unsigned int _nKernelParams;

    // Set the Kernel's ideas about the x- and y- coordinates
    virtual void _setKernelXY() {}
};

typedef std::vector<std::shared_ptr<Kernel>> KernelList;

class FixedKernel : public afw::table::io::PersistableFacade<FixedKernel>, public Kernel {
public:
    explicit FixedKernel();

    explicit FixedKernel(lsst::afw::image::Image<Pixel> const &image  
    );

    explicit FixedKernel(lsst::afw::math::Kernel const &kernel,  
                         lsst::geom::Point2D const &pos          
    );

    FixedKernel(const FixedKernel &) = delete;
    FixedKernel(FixedKernel &&) = delete;
    FixedKernel &operator=(const FixedKernel &) = delete;
    FixedKernel &operator=(FixedKernel &&) = delete;

    ~FixedKernel() override = default;

    std::shared_ptr<Kernel> clone() const override;

    std::shared_ptr<Kernel> resized(int width, int height) const override;

    std::string toString(std::string const &prefix = "") const override;

    virtual Pixel getSum() const { return _sum; }

    bool isPersistable() const noexcept override { return true; }

    class Factory;

protected:
    double doComputeImage(lsst::afw::image::Image<Pixel> &image, bool doNormalize) const override;

    std::string getPersistenceName() const override;

    void write(OutputArchiveHandle &handle) const override;

private:
    lsst::afw::image::Image<Pixel> _image;
    Pixel _sum;
};

class AnalyticKernel : public afw::table::io::PersistableFacade<AnalyticKernel>, public Kernel {
public:
    typedef lsst::afw::math::Function2<Pixel> KernelFunction;
    typedef std::shared_ptr<lsst::afw::math::Function2<Pixel>> KernelFunctionPtr;

    explicit AnalyticKernel();

    explicit AnalyticKernel(int width, int height, KernelFunction const &kernelFunction,
                            Kernel::SpatialFunction const &spatialFunction = NullSpatialFunction());

    explicit AnalyticKernel(int width, int height, KernelFunction const &kernelFunction,
                            std::vector<Kernel::SpatialFunctionPtr> const &spatialFunctionList);

    AnalyticKernel(const AnalyticKernel &) = delete;
    AnalyticKernel(AnalyticKernel &&) = delete;
    AnalyticKernel &operator=(const AnalyticKernel &) = delete;
    AnalyticKernel &operator=(AnalyticKernel &&) = delete;

    ~AnalyticKernel() override = default;

    std::shared_ptr<Kernel> clone() const override;

    std::shared_ptr<Kernel> resized(int width, int height) const override;

    double computeImage(lsst::afw::image::Image<Pixel> &image, bool doNormalize, double x = 0.0,
                        double y = 0.0) const;

    std::vector<double> getKernelParameters() const override;

    virtual KernelFunctionPtr getKernelFunction() const;

    std::string toString(std::string const &prefix = "") const override;

    bool isPersistable() const noexcept override { return true; }

    class Factory;

protected:
    double doComputeImage(lsst::afw::image::Image<Pixel> &image, bool doNormalize) const override;

    std::string getPersistenceName() const override;

    void write(OutputArchiveHandle &handle) const override;

protected:
    void setKernelParameter(unsigned int ind, double value) const override;

    KernelFunctionPtr _kernelFunctionPtr;
};

class DeltaFunctionKernel : public afw::table::io::PersistableFacade<DeltaFunctionKernel>, public Kernel {
public:
    // Traits values for this class of Kernel
    typedef deltafunction_kernel_tag kernel_fill_factor;

    explicit DeltaFunctionKernel(int width, int height, lsst::geom::Point2I const &point);

    DeltaFunctionKernel(const DeltaFunctionKernel &) = delete;
    DeltaFunctionKernel(DeltaFunctionKernel &&) = delete;
    DeltaFunctionKernel &operator=(const DeltaFunctionKernel &) = delete;
    DeltaFunctionKernel &operator=(DeltaFunctionKernel &&) = delete;

    ~DeltaFunctionKernel() override = default;

    std::shared_ptr<Kernel> clone() const override;

    std::shared_ptr<Kernel> resized(int width, int height) const override;

    lsst::geom::Point2I getPixel() const { return _pixel; }

    std::string toString(std::string const &prefix = "") const override;

    bool isPersistable() const noexcept override { return true; }

    class Factory;

protected:
    double doComputeImage(lsst::afw::image::Image<Pixel> &image, bool doNormalize) const override;

    std::string getPersistenceName() const override;

    void write(OutputArchiveHandle &handle) const override;

private:
    lsst::geom::Point2I _pixel;
};

class LinearCombinationKernel : public afw::table::io::PersistableFacade<LinearCombinationKernel>,
                                public Kernel {
public:
    explicit LinearCombinationKernel();

    explicit LinearCombinationKernel(KernelList const &kernelList,
                                     std::vector<double> const &kernelParameters);

    explicit LinearCombinationKernel(KernelList const &kernelList,
                                     Kernel::SpatialFunction const &spatialFunction);

    explicit LinearCombinationKernel(KernelList const &kernelList,
                                     std::vector<Kernel::SpatialFunctionPtr> const &spatialFunctionList);

    LinearCombinationKernel(const LinearCombinationKernel &) = delete;
    LinearCombinationKernel(LinearCombinationKernel &&) = delete;
    LinearCombinationKernel &operator=(const LinearCombinationKernel &) = delete;
    LinearCombinationKernel &operator=(LinearCombinationKernel &&) = delete;

    ~LinearCombinationKernel() override = default;

    std::shared_ptr<Kernel> clone() const override;

    std::shared_ptr<Kernel> resized(int width, int height) const override;

    std::vector<double> getKernelParameters() const override;

    virtual KernelList const &getKernelList() const;

    std::vector<double> getKernelSumList() const;

    int getNBasisKernels() const { return static_cast<int>(_kernelList.size()); };

    void checkKernelList(const KernelList &kernelList) const;

    bool isDeltaFunctionBasis() const { return _isDeltaFunctionBasis; };

    std::shared_ptr<Kernel> refactor() const;

    std::string toString(std::string const &prefix = "") const override;

    bool isPersistable() const noexcept override { return true; }

    class Factory;

protected:
    double doComputeImage(lsst::afw::image::Image<Pixel> &image, bool doNormalize) const override;

    std::string getPersistenceName() const override;

    void write(OutputArchiveHandle &handle) const override;

    void setKernelParameter(unsigned int ind, double value) const override;

private:
    void _setKernelList(KernelList const &kernelList);

    KernelList _kernelList;  
    std::vector<std::shared_ptr<lsst::afw::image::Image<Pixel>>> _kernelImagePtrList;
    std::vector<double> _kernelSumList;  
    mutable std::vector<double> _kernelParams;
    bool _isDeltaFunctionBasis;
};

class SeparableKernel : public afw::table::io::PersistableFacade<SeparableKernel>, public Kernel {
public:
    typedef lsst::afw::math::Function1<Pixel> KernelFunction;
    typedef std::shared_ptr<KernelFunction> KernelFunctionPtr;

    explicit SeparableKernel();

    explicit SeparableKernel(int width, int height, KernelFunction const &kernelColFunction,
                             KernelFunction const &kernelRowFunction,
                             Kernel::SpatialFunction const &spatialFunction = NullSpatialFunction());

    explicit SeparableKernel(int width, int height, KernelFunction const &kernelColFunction,
                             KernelFunction const &kernelRowFunction,
                             std::vector<Kernel::SpatialFunctionPtr> const &spatialFunctionList);

    SeparableKernel(const SeparableKernel &) = delete;
    SeparableKernel(SeparableKernel &&) = delete;
    SeparableKernel &operator=(const SeparableKernel &) = delete;
    SeparableKernel &operator=(SeparableKernel &&) = delete;

    ~SeparableKernel() override = default;

    std::shared_ptr<Kernel> clone() const override;

    std::shared_ptr<Kernel> resized(int width, int height) const override;

    double computeVectors(std::vector<Pixel> &colList, std::vector<Pixel> &rowList, bool doNormalize,
                          double x = 0.0, double y = 0.0) const;

    double getKernelParameter(unsigned int i) const override {
        unsigned int const ncol = _kernelColFunctionPtr->getNParameters();
        if (i < ncol) {
            return _kernelColFunctionPtr->getParameter(i);
        } else {
            i -= ncol;
            return _kernelRowFunctionPtr->getParameter(i);
        }
    }
    std::vector<double> getKernelParameters() const override;

    KernelFunctionPtr getKernelColFunction() const;

    KernelFunctionPtr getKernelRowFunction() const;

    std::string toString(std::string const &prefix = "") const override;

    /***
     * Compute a cache of values for the x and y kernel functions
     *
     * A value of 0 disables the cache for maximum accuracy.
     * 10,000 typically results in a warping error of a fraction of a count.
     * 100,000 typically results in a warping error of less than 0.01 count.
     *
     * @param cacheSize cache size (number of double precision array elements in the x and y caches)
     */
    void computeCache(int const cacheSize) override;

    int getCacheSize() const override;

protected:
    double doComputeImage(lsst::afw::image::Image<Pixel> &image, bool doNormalize) const override;

    void setKernelParameter(unsigned int ind, double value) const override;

private:
    double basicComputeVectors(std::vector<Pixel> &colList, std::vector<Pixel> &rowList,
                               bool doNormalize) const;

    KernelFunctionPtr _kernelColFunctionPtr;
    KernelFunctionPtr _kernelRowFunctionPtr;
    mutable std::vector<Pixel> _localColList;  // used by doComputeImage
    mutable std::vector<Pixel> _localRowList;
    mutable std::vector<double> _kernelX;  // used by SeparableKernel::basicComputeVectors
    mutable std::vector<double> _kernelY;
    //
    // Cached values of the row- and column- kernels
    //
    mutable std::vector<std::vector<double>> _kernelRowCache;
    mutable std::vector<std::vector<double>> _kernelColCache;

    virtual void _setKernelXY() override {
        lsst::geom::Extent2I const dim = getDimensions();
        lsst::geom::Point2I const ctr = getCtr();

        assert(dim[0] == static_cast<int>(_kernelX.size()));
        for (int i = 0; i != dim.getX(); ++i) {
            _kernelX[i] = i - ctr.getX();
        }

        assert(dim[1] == static_cast<int>(_kernelY.size()));
        for (int i = 0; i != dim.getY(); ++i) {
            _kernelY[i] = i - ctr.getY();
        }
    }
};
}  // namespace math
}  // namespace afw
}  // namespace lsst

#endif  // !defined(LSST_AFW_MATH_KERNEL_H)