Convolve.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_DETAIL_CONVOLVE_H
#define LSST_AFW_MATH_DETAIL_CONVOLVE_H

#include <memory>
#include <sstream>

#include "lsst/pex/exceptions.h"
#include "lsst/afw/geom.h"
#include "lsst/afw/image/Image.h"
#include "lsst/afw/image/MaskedImage.h"
#include "lsst/afw/math/Kernel.h"
#include "lsst/afw/math/ConvolveImage.h"

#define IS_INSTANCE(A, B) (dynamic_cast<B const*>(&(A)) != NULL)

namespace lsst {
namespace afw {
namespace math {

namespace detail {
    template <typename OutImageT, typename InImageT>
    void basicConvolve(
            OutImageT& convolvedImage,
            InImageT const& inImage,
            lsst::afw::math::Kernel const& kernel,
            lsst::afw::math::ConvolutionControl const& convolutionControl);

    template <typename OutImageT, typename InImageT>
    void basicConvolve(
            OutImageT& convolvedImage,
            InImageT const& inImage,
            lsst::afw::math::DeltaFunctionKernel const& kernel,
            lsst::afw::math::ConvolutionControl const&);

    template <typename OutImageT, typename InImageT>
    void basicConvolve(
            OutImageT& convolvedImage,
            InImageT const& inImage,
            lsst::afw::math::LinearCombinationKernel const& kernel,
            lsst::afw::math::ConvolutionControl const& convolutionControl);

    template <typename OutImageT, typename InImageT>
    void basicConvolve(
            OutImageT& convolvedImage,
            InImageT const& inImage,
            lsst::afw::math::SeparableKernel const& kernel,
            lsst::afw::math::ConvolutionControl const& convolutionControl);

    template <typename OutImageT, typename InImageT>
    void convolveWithBruteForce(
            OutImageT &convolvedImage,
            InImageT const& inImage,
            lsst::afw::math::Kernel const& kernel,
            lsst::afw::math::ConvolutionControl const& convolutionControl);

    // I would prefer this to be nested in KernelImagesForRegion but SWIG doesn't support that
    class RowOfKernelImagesForRegion;

    class KernelImagesForRegion :
        public lsst::daf::base::Citizen,
        public lsst::daf::base::Persistable
    {
    public:
        typedef CONST_PTR(lsst::afw::math::Kernel) KernelConstPtr;
        typedef lsst::afw::image::Image<lsst::afw::math::Kernel::Pixel> Image;
        typedef PTR(Image) ImagePtr;
        typedef CONST_PTR(Image) ImageConstPtr;
        typedef CONST_PTR(KernelImagesForRegion) ConstPtr;
        typedef PTR(KernelImagesForRegion) Ptr;

        enum Location {
            BOTTOM_LEFT, BOTTOM_RIGHT, TOP_LEFT, TOP_RIGHT
        };

        KernelImagesForRegion(
                KernelConstPtr kernelPtr,
                lsst::afw::geom::Box2I const &bbox,
                lsst::afw::geom::Point2I const &xy0,
                bool doNormalize);
        KernelImagesForRegion(
                KernelConstPtr kernelPtr,
                lsst::afw::geom::Box2I const &bbox,
                lsst::afw::geom::Point2I const &xy0,
                bool doNormalize,
                ImagePtr bottomLeftImagePtr,
                ImagePtr bottomRightImagePtr,
                ImagePtr topLeftImagePtr,
                ImagePtr topRightImagePtr);

        lsst::afw::geom::Box2I getBBox() const { return _bbox; };
        lsst::afw::geom::Point2I getXY0() const { return _xy0; };
        bool getDoNormalize() const { return _doNormalize; };
        ImagePtr getImage(Location location) const;
        KernelConstPtr getKernel() const { return _kernelPtr; };
        lsst::afw::geom::Point2I getPixelIndex(Location location) const;
        bool computeNextRow(RowOfKernelImagesForRegion &regionRow) const;

        static int getMinInterpolationSize() { return _MinInterpolationSize; };
    private:
        typedef std::vector<Location> LocationList;

        void _computeImage(Location location) const;
        inline void _insertImage(Location location, ImagePtr imagePtr) const;
        void _moveUp(bool isFirst, int newHeight);

        // static helper functions
        static inline int _computeNextSubregionLength(int length, int nDivisions);
        static std::vector<int> _computeSubregionLengths(int length, int nDivisions);

        // member variables
        KernelConstPtr _kernelPtr;
        lsst::afw::geom::Box2I _bbox;
        lsst::afw::geom::Point2I _xy0;
        bool _doNormalize;
        mutable std::vector<ImagePtr> _imagePtrList;

        static int const _MinInterpolationSize;
    };

    class RowOfKernelImagesForRegion {
    public:
        typedef std::vector<PTR(KernelImagesForRegion)> RegionList;
        typedef RegionList::iterator Iterator;
        typedef RegionList::const_iterator ConstIterator;

        RowOfKernelImagesForRegion(int nx, int ny);
        RegionList::const_iterator begin() const { return _regionList.begin(); };
        RegionList::const_iterator end() const { return _regionList.end(); };
        RegionList::iterator begin() { return _regionList.begin(); };
        RegionList::iterator end() { return _regionList.end(); };
        PTR(KernelImagesForRegion) front() { return _regionList.front(); };
        PTR(KernelImagesForRegion) back() { return _regionList.back(); };
        int getNX() const { return _nx; };
        int getNY() const { return _ny; };
        int getYInd() const { return _yInd; };
        CONST_PTR(KernelImagesForRegion) getRegion(int ind) const { return _regionList.at(ind); };
        bool hasData() const { return static_cast<bool>(_regionList[0]); };
        bool isLastRow() const { return _yInd + 1 >= _ny; };
        int incrYInd() { return ++_yInd; };

    private:
        int _nx;
        int _ny;
        int _yInd;
        RegionList _regionList;
    };

    template <typename OutImageT, typename InImageT>
    void convolveWithInterpolation(
            OutImageT &outImage,
            InImageT const &inImage,
            lsst::afw::math::Kernel const &kernel,
            ConvolutionControl const &convolutionControl);

    struct ConvolveWithInterpolationWorkingImages {
    public:
        typedef lsst::afw::image::Image<lsst::afw::math::Kernel::Pixel> Image;
        ConvolveWithInterpolationWorkingImages(geom::Extent2I const & dimensions) :
            leftImage(dimensions),
            rightImage(dimensions),
            leftDeltaImage(dimensions),
            rightDeltaImage(dimensions),
            deltaImage(dimensions),
            kernelImage(dimensions)
        { }
        Image leftImage;
        Image rightImage;
        Image leftDeltaImage;
        Image rightDeltaImage;
        Image deltaImage;
        Image kernelImage;
    };

    template <typename OutImageT, typename InImageT>
    void convolveRegionWithInterpolation(
            OutImageT &outImage,
            InImageT const &inImage,
            KernelImagesForRegion const &region,
            ConvolveWithInterpolationWorkingImages &workingImages);
}}}}   // lsst::afw::math::detail

/*
 * Define inline functions
 */

inline int lsst::afw::math::detail::KernelImagesForRegion::_computeNextSubregionLength(
    int length,     
    int nDivisions) 
{
    return static_cast<int>(std::floor(0.5 +
        (static_cast<double>(length) / static_cast<double>(nDivisions))));
}

inline void lsst::afw::math::detail::KernelImagesForRegion::_insertImage(
        Location location,      
        ImagePtr imagePtr)      
const {
    if (imagePtr) {
        if (_kernelPtr->getDimensions() != imagePtr->getDimensions()) {
            std::ostringstream os;
            os << "image dimensions = ( "
                << imagePtr->getWidth() << ", " << imagePtr->getHeight()
                << ") != (" << _kernelPtr->getWidth() << ", " << _kernelPtr->getHeight()
                << ") = kernel dimensions";
            throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterError, os.str());
        }
        _imagePtrList[location] = imagePtr;
    }
}

#endif // !defined(LSST_AFW_MATH_DETAIL_CONVOLVE_H)