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/format.hpp"
 
#include "lsst/geom.h"
#include "lsst/afw/image/Image.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 afw::table::io::PersistableFacade<Kernel>, public afw::table::io::Persistable {
public:
    using Pixel = double;
    using SpatialFunctionPtr = std::shared_ptr<lsst::afw::math::Function2<double>>;
    using SpatialFunction = lsst::afw::math::Function2<double>;
    using NullSpatialFunction = lsst::afw::math::NullFunction2<double>;
 
    // Traits values for this class of Kernel
    using kernel_fill_factor = generic_kernel_tag;
 
    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 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 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() {}
};
 
using KernelList = std::vector<std::shared_ptr<Kernel>>;
 
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:
    using KernelFunction = lsst::afw::math::Function2<Pixel>;
    using KernelFunctionPtr = std::shared_ptr<lsst::afw::math::Function2<Pixel>>;
 
    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
    using kernel_fill_factor = deltafunction_kernel_tag;
 
    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:
    using KernelFunction = lsst::afw::math::Function1<Pixel>;
    using KernelFunctionPtr = std::shared_ptr<KernelFunction>;
 
    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)