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
/*
 * Convolution support
 */
#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& convolutionControl);

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:
    typedef std::shared_ptr<lsst::afw::math::Kernel const> KernelConstPtr;
    typedef lsst::afw::image::Image<lsst::afw::math::Kernel::Pixel> Image;
    typedef std::shared_ptr<Image> ImagePtr;
    typedef std::shared_ptr<Image const> ImageConstPtr;

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

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

    lsst::geom::Box2I getBBox() const { return _bbox; };
    lsst::geom::Point2I getXY0() const { return _xy0; };
    bool getDoNormalize() const { return _doNormalize; };
    ImagePtr getImage(Location location) const;
    KernelConstPtr getKernel() const { return _kernelPtr; };
    lsst::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::geom::Box2I _bbox;
    lsst::geom::Point2I _xy0;
    bool _doNormalize;
    mutable std::vector<ImagePtr> _imagePtrList;

    static int const _MinInterpolationSize;
};

class RowOfKernelImagesForRegion final {
public:
    typedef std::vector<std::shared_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(); };
    std::shared_ptr<KernelImagesForRegion> front() { return _regionList.front(); };
    std::shared_ptr<KernelImagesForRegion> back() { return _regionList.back(); };
    int getNX() const { return _nx; };
    int getNY() const { return _ny; };
    int getYInd() const { return _yInd; };
    std::shared_ptr<KernelImagesForRegion const> 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 final {
public:
    typedef lsst::afw::image::Image<lsst::afw::math::Kernel::Pixel> Image;
    ConvolveWithInterpolationWorkingImages(lsst::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);

/*
 * Define inline functions
 */

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

inline void 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;
    }
}
}  // namespace detail
}  // namespace math
}  // namespace afw
}  // namespace lsst

#endif  // !defined(LSST_AFW_MATH_DETAIL_CONVOLVE_H)