Kernel.h

Go to the documentation of this file.

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

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

#include <utility>
#include <vector>

#include "boost/mpl/or.hpp"
#include "boost/shared_ptr.hpp"
#include "boost/make_shared.hpp"
#include "boost/static_assert.hpp"
#include "boost/type_traits/is_same.hpp"
#include "boost/type_traits/is_base_and_derived.hpp"

#include "boost/serialization/shared_ptr.hpp"
#include "boost/serialization/vector.hpp"
#include "boost/serialization/export.hpp"

#include "lsst/daf/base/Persistable.h"
#include "lsst/daf/base/Citizen.h"
#include "lsst/afw/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 formatters {
class KernelFormatter;
}

namespace math {

#ifndef SWIG
using boost::serialization::make_nvp;
#endif

    class Kernel : public lsst::daf::base::Citizen, public lsst::daf::base::Persistable,
                   public afw::table::io::PersistableFacade<Kernel>,
                   public afw::table::io::Persistable
    {

    public:
        typedef double Pixel;
        typedef PTR(Kernel) Ptr;
        typedef CONST_PTR(Kernel) ConstPtr;
        typedef 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
        );

        virtual ~Kernel() {}

        virtual PTR(Kernel) clone() const = 0;

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

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

        void setDimensions(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::afw::geom::Point2I getCtr() const {
            return lsst::afw::geom::Point2I(_ctrX, _ctrY);
        }

        inline int getCtrX() const {
            return _ctrX;
        }

        inline int getCtrY() const {
            return _ctrY;
        }
        
        inline lsst::afw::geom::Box2I getBBox() const {
            return lsst::afw::geom::Box2I(
                lsst::afw::geom::Point2I(-_ctrX, -_ctrY),
                lsst::afw::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::afw::geom::Box2I growBBox(lsst::afw::geom::Box2I const &bbox) const;
        
        lsst::afw::geom::Box2I shrinkBBox(lsst::afw::geom::Box2I const &bbox) const;

        inline void setCtr(lsst::afw::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:

        virtual std::string getPythonModule() const;

        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:
        LSST_PERSIST_FORMATTER(lsst::afw::formatters::KernelFormatter)

        int _width;
        int _height;
        int _ctrX;
        int _ctrY;
        unsigned int _nKernelParams;
        
        // prevent copying and assignment (to avoid problems from type slicing)
        Kernel(const Kernel&);
        Kernel& operator=(const Kernel&);
        // Set the Kernel's ideas about the x- and y- coordinates
        virtual void _setKernelXY() {}
    };

    typedef std::vector<PTR(Kernel)> KernelList;

    class FixedKernel : public afw::table::io::PersistableFacade<FixedKernel>, public Kernel {
    public:
        typedef PTR(FixedKernel) Ptr;
        typedef CONST_PTR(FixedKernel) ConstPtr;

        explicit FixedKernel();

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

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

        virtual ~FixedKernel() {}

        virtual PTR(Kernel) clone() const;

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

        virtual Pixel getSum() const {
            return _sum;
        }

        virtual bool isPersistable() const { return true; }

        class Factory;

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

        virtual std::string getPersistenceName() const;

        virtual void write(OutputArchiveHandle & handle) const;

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

        friend class boost::serialization::access;
        template <class Archive>
            void serialize(Archive& ar, unsigned int const version) {
                ar & make_nvp("k",
                        boost::serialization::base_object<Kernel>(*this));
                ar & make_nvp("img", _image);
                ar & make_nvp("sum", _sum);
            }
    };


    class AnalyticKernel : public afw::table::io::PersistableFacade<AnalyticKernel>, public Kernel {
    public:
        typedef PTR(AnalyticKernel) Ptr;
        typedef CONST_PTR(AnalyticKernel) ConstPtr;
        typedef lsst::afw::math::Function2<Pixel> KernelFunction;
        typedef 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 
        );

        virtual ~AnalyticKernel() {}

        virtual PTR(Kernel) clone() const;

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

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

        virtual KernelFunctionPtr getKernelFunction() const;

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

        virtual bool isPersistable() const { return true; }

        class Factory;

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

        virtual std::string getPersistenceName() const;

        virtual void write(OutputArchiveHandle & handle) const;

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

        KernelFunctionPtr _kernelFunctionPtr;

        friend class boost::serialization::access;
        template <class Archive>
            void serialize(Archive& ar, unsigned int const version) {
                ar & make_nvp("k",
                        boost::serialization::base_object<Kernel>(*this));
                ar & make_nvp("fn", _kernelFunctionPtr);
            }
    };


    class DeltaFunctionKernel : public afw::table::io::PersistableFacade<DeltaFunctionKernel>,
                                public Kernel
    {
    public:
        typedef PTR(DeltaFunctionKernel) Ptr;
        typedef CONST_PTR(DeltaFunctionKernel) ConstPtr;
        // Traits values for this class of Kernel
        typedef deltafunction_kernel_tag kernel_fill_factor;

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

        virtual ~DeltaFunctionKernel() {}

        virtual PTR(Kernel) clone() const;

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

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

        virtual bool isPersistable() const { return true; }

        class Factory;

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

        virtual std::string getPersistenceName() const;

        virtual void write(OutputArchiveHandle & handle) const;

    private:
        lsst::afw::geom::Point2I _pixel;

        friend class boost::serialization::access;
        template <class Archive>
        void serialize(Archive& ar, unsigned int const version) {
            boost::serialization::void_cast_register<
                DeltaFunctionKernel, Kernel>(
                    static_cast<DeltaFunctionKernel*>(0),
                    static_cast<Kernel*>(0));
        }
    };


    class LinearCombinationKernel : public afw::table::io::PersistableFacade<LinearCombinationKernel>,
                                    public Kernel
    {
    public:
        typedef PTR(LinearCombinationKernel) Ptr;
        typedef CONST_PTR(LinearCombinationKernel) ConstPtr;

        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
        );

        virtual ~LinearCombinationKernel() {}

        virtual PTR(Kernel) clone() const;

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

        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; };
        
        PTR(Kernel) refactor() const;

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

        virtual bool isPersistable() const { return true; }

        class Factory;

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

        virtual std::string getPersistenceName() const;

        virtual void write(OutputArchiveHandle & handle) const;

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

    private:
        void _setKernelList(KernelList const &kernelList);
        
        KernelList _kernelList; 
        std::vector<PTR(lsst::afw::image::Image<Pixel>)> _kernelImagePtrList;
        std::vector<double> _kernelSumList; 
        mutable std::vector<double> _kernelParams;
        bool _isDeltaFunctionBasis;

        friend class boost::serialization::access;
        template <class Archive>
            void serialize(Archive& ar, unsigned int const version) {
                ar & make_nvp("k", boost::serialization::base_object<Kernel>(*this));
                ar & make_nvp("klist", _kernelList);
                ar & make_nvp("kimglist", _kernelImagePtrList);
                ar & make_nvp("ksumlist", _kernelSumList);
                ar & make_nvp("params", _kernelParams);
                if (version > 0) {
                    ar & make_nvp("deltaBasis", _isDeltaFunctionBasis);
                }
                else if (Archive::is_loading::value) {
                    _isDeltaFunctionBasis = false;
                }
            }
    };

    class SeparableKernel : public afw::table::io::PersistableFacade<SeparableKernel>, public Kernel {
    public:
        typedef PTR(SeparableKernel) Ptr;
        typedef CONST_PTR(SeparableKernel) ConstPtr;
        typedef lsst::afw::math::Function1<Pixel> KernelFunction;
        typedef 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 
        );
        virtual ~SeparableKernel() {}

        virtual PTR(Kernel) clone() const;

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

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

        KernelFunctionPtr getKernelColFunction() const;

        KernelFunctionPtr getKernelRowFunction() const;

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

        /***
         * @brief 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.
         */
        virtual void computeCache(
            int const cacheSize 
        );
        
        virtual int getCacheSize() const;

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

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

    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;

        friend class boost::serialization::access;
        template <class Archive>
            void serialize(Archive& ar, unsigned int const version) {
                ar & make_nvp("k",
                    boost::serialization::base_object<Kernel>(*this));
                ar & make_nvp("colfn", _kernelColFunctionPtr);
                ar & make_nvp("rowfn", _kernelRowFunctionPtr);
                ar & make_nvp("cols", _localColList);
                ar & make_nvp("rows", _localRowList);
                ar & make_nvp("kernelX", _kernelX);
                ar & make_nvp("kernelY", _kernelY);
            }

        virtual void _setKernelXY() {
            lsst::afw::geom::Extent2I const dim = getDimensions();
            lsst::afw::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();
            }
        }
    };

}}}   // lsst:afw::math

namespace boost {
namespace serialization {

template <class Archive>
inline void save_construct_data(
    Archive& ar, lsst::afw::math::DeltaFunctionKernel const* k,
    unsigned int const file_version) {
    int width = k->getWidth();
    int height = k->getHeight();
    int x = k->getPixel().getX();
    int y = k->getPixel().getY();
    ar << make_nvp("width", width);
    ar << make_nvp("height", height);
    ar << make_nvp("pixX", x);
    ar << make_nvp("pixY", y);
}

template <class Archive>
inline void load_construct_data(
    Archive& ar, lsst::afw::math::DeltaFunctionKernel* k,
    unsigned int const file_version) {
    int width;
    int height;
    int x;
    int y;
    ar >> make_nvp("width", width);
    ar >> make_nvp("height", height);
    ar >> make_nvp("pixX", x);
    ar >> make_nvp("pixY", y);
    ::new(k) lsst::afw::math::DeltaFunctionKernel(
        width, height, lsst::afw::geom::Point2I(x, y));
}

}}

#ifndef SWIG
BOOST_CLASS_VERSION(lsst::afw::math::LinearCombinationKernel, 1)
#endif

#endif // !defined(LSST_AFW_MATH_KERNEL_H)