ConvolveImage.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_CONVOLVEIMAGE_H
#define LSST_AFW_MATH_CONVOLVEIMAGE_H
/*
 * Convolve and convolveAtAPoint functions for Image and Kernel
 *
 * @todo Consider adding a flag to convolve indicating which specialized version of basicConvolve was used.
 *   This would only be used for unit testing and trace messages suffice (barely), so not a high priority.
 */
#include <limits>
#include <sstream>

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

namespace lsst {
namespace afw {
namespace math {

class ConvolutionControl {
public:
    ConvolutionControl(bool doNormalize = true,  
                       bool doCopyEdge = false,  
                       int maxInterpolationDistance = 10  
                       )
            : _doNormalize(doNormalize),
              _doCopyEdge(doCopyEdge),
              _maxInterpolationDistance(maxInterpolationDistance) {}

    bool getDoNormalize() const { return _doNormalize; }
    bool getDoCopyEdge() const { return _doCopyEdge; }
    int getMaxInterpolationDistance() const { return _maxInterpolationDistance; };

    void setDoNormalize(bool doNormalize) { _doNormalize = doNormalize; }
    void setDoCopyEdge(bool doCopyEdge) { _doCopyEdge = doCopyEdge; }
    void setMaxInterpolationDistance(int maxInterpolationDistance) {
        _maxInterpolationDistance = maxInterpolationDistance;
    }

private:
    bool _doNormalize;              
    bool _doCopyEdge;               
    int _maxInterpolationDistance;  
};

template <typename OutImageT, typename InImageT>
void scaledPlus(OutImageT& outImage, double c1, InImageT const& inImage1, double c2,
                InImageT const& inImage2);

template <typename OutImageT, typename InImageT>
inline typename OutImageT::SinglePixel convolveAtAPoint(
        typename InImageT::const_xy_locator inImageLocator,
        typename lsst::afw::image::Image<lsst::afw::math::Kernel::Pixel>::const_xy_locator kernelLocator,
        int kWidth, int kHeight);

template <typename OutImageT, typename InImageT>
inline typename OutImageT::SinglePixel convolveAtAPoint(
        typename InImageT::const_xy_locator inImageLocator,
        std::vector<lsst::afw::math::Kernel::Pixel> const& kernelColList,
        std::vector<lsst::afw::math::Kernel::Pixel> const& kernelRowList);

template <typename OutImageT, typename InImageT, typename KernelT>
void convolve(OutImageT& convolvedImage, InImageT const& inImage, KernelT const& kernel,
              ConvolutionControl const& convolutionControl = ConvolutionControl());

template <typename OutImageT, typename InImageT, typename KernelT>
void convolve(OutImageT& convolvedImage, InImageT const& inImage, KernelT const& kernel, bool doNormalize,
              bool doCopyEdge = false);

template <typename ImageT>
typename ImageT::SinglePixel edgePixel(lsst::afw::image::detail::Image_tag
) {
    typedef typename ImageT::SinglePixel SinglePixelT;
    return SinglePixelT(std::numeric_limits<SinglePixelT>::has_quiet_NaN
                                ? std::numeric_limits<SinglePixelT>::quiet_NaN()
                                : 0);
}

template <typename MaskedImageT>
typename MaskedImageT::SinglePixel edgePixel(
        lsst::afw::image::detail::MaskedImage_tag
) {
    typedef typename MaskedImageT::Image::Pixel ImagePixelT;
    typedef typename MaskedImageT::Variance::Pixel VariancePixelT;

    auto imagePixel = std::numeric_limits<ImagePixelT>::has_quiet_NaN
                              ? std::numeric_limits<ImagePixelT>::quiet_NaN()
                              : 0;
    auto maskPixel = MaskedImageT::Mask::getPlaneBitMask("NO_DATA");
    auto variancePixel = std::numeric_limits<VariancePixelT>::has_infinity
                                 ? std::numeric_limits<VariancePixelT>::infinity()
                                 : 0;
    return typename MaskedImageT::SinglePixel(imagePixel, maskPixel, variancePixel);
}

/*
 * Define inline functions
 */

template <typename OutImageT, typename InImageT>
inline typename OutImageT::SinglePixel convolveAtAPoint(
        typename InImageT::const_xy_locator inImageLocator,
        image::Image<Kernel::Pixel>::const_xy_locator kernelLocator, int kWidth, int kHeight) {
    typename OutImageT::SinglePixel outValue = 0;
    for (int kRow = 0; kRow != kHeight; ++kRow) {
        for (image::Image<Kernel::Pixel>::const_xy_locator
                     kEnd = kernelLocator + image::detail::difference_type(kWidth, 0);
             kernelLocator != kEnd; ++inImageLocator.x(), ++kernelLocator.x()) {
            typename Kernel::Pixel const kVal = kernelLocator[0];
            if (kVal != 0) {
                outValue += *inImageLocator * kVal;
            }
        }

        inImageLocator += image::detail::difference_type(-kWidth, 1);
        kernelLocator += image::detail::difference_type(-kWidth, 1);
    }

    return outValue;
}

template <typename OutImageT, typename InImageT>
inline typename OutImageT::SinglePixel convolveAtAPoint(typename InImageT::const_xy_locator inImageLocator,
                                                        std::vector<Kernel::Pixel> const& kernelXList,
                                                        std::vector<Kernel::Pixel> const& kernelYList) {
    typedef typename std::vector<Kernel::Pixel>::const_iterator k_iter;

    typedef typename OutImageT::SinglePixel OutT;
    OutT outValue = 0;
    for (k_iter kernelYIter = kernelYList.begin(), yEnd = kernelYList.end(); kernelYIter != yEnd;
         ++kernelYIter) {
        OutT outValueY = 0;
        for (k_iter kernelXIter = kernelXList.begin(), xEnd = kernelXList.end(); kernelXIter != xEnd;
             ++kernelXIter, ++inImageLocator.x()) {
            typename Kernel::Pixel const kValX = *kernelXIter;
            if (kValX != 0) {
                outValueY += *inImageLocator * kValX;
            }
        }

        double const kValY = *kernelYIter;
        if (kValY != 0) {
            outValue += outValueY * kValY;
        }

        inImageLocator += image::detail::difference_type(-kernelXList.size(), 1);
    }

    return outValue;
}
}  // namespace math
}  // namespace afw
}  // namespace lsst

#endif  // !defined(LSST_AFW_MATH_CONVOLVEIMAGE_H)