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:
    using KernelConstPtr = std::shared_ptr<const lsst::afw::math::Kernel>;
    using Image = lsst::afw::image::Image<lsst::afw::math::Kernel::Pixel>;
    using ImagePtr = std::shared_ptr<Image>;
    using ImageConstPtr = std::shared_ptr<const Image>;
 
    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:
    using LocationList = std::vector<Location>;
 
    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:
    using RegionList = std::vector<std::shared_ptr<KernelImagesForRegion>>;
    using Iterator = RegionList::iterator;
    using ConstIterator = RegionList::const_iterator;
 
    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:
    using Image = lsst::afw::image::Image<lsst::afw::math::Kernel::Pixel>;
    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)