WarpAtOnePoint.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/>.
 */
#include <vector>
 
#include "lsst/afw/math/Kernel.h"
#include "lsst/afw/image/Image.h"
#include "lsst/afw/image/MaskedImage.h"
#include "lsst/geom/Point.h"
 
namespace lsst {
namespace afw {
namespace math {
namespace detail {
 
template <typename DestImageT, typename SrcImageT>
class WarpAtOnePoint final {
public:
    WarpAtOnePoint(SrcImageT const &srcImage, WarpingControl const &control,
                   typename DestImageT::SinglePixel padValue)
            : _srcImage(srcImage),
              _kernelPtr(control.getWarpingKernel()),
              _maskKernelPtr(control.getMaskWarpingKernel()),
              _hasMaskKernel(control.getMaskWarpingKernel()),
              _kernelCtr(_kernelPtr->getCtr()),
              _maskKernelCtr(_maskKernelPtr ? _maskKernelPtr->getCtr() : lsst::geom::Point2I(0, 0)),
              _growFullMask(control.getGrowFullMask()),
              _xList(_kernelPtr->getWidth()),
              _yList(_kernelPtr->getHeight()),
              _maskXList(_maskKernelPtr ? _maskKernelPtr->getWidth() : 0),
              _maskYList(_maskKernelPtr ? _maskKernelPtr->getHeight() : 0),
              _padValue(padValue),
              _srcGoodBBox(_kernelPtr->shrinkBBox(srcImage.getBBox(lsst::afw::image::LOCAL))){};
 
    bool operator()(typename DestImageT::x_iterator &destXIter, lsst::geom::Point2D const &srcPos,
                    double relativeArea, lsst::afw::image::detail::Image_tag) {
        // Compute associated source pixel index as integer and nonnegative fractional parts;
        // the latter is used to compute the remapping kernel.
        std::pair<int, double> srcIndFracX = _srcImage.positionToIndex(srcPos[0], lsst::afw::image::X);
        std::pair<int, double> srcIndFracY = _srcImage.positionToIndex(srcPos[1], lsst::afw::image::Y);
        if (srcIndFracX.second < 0) {
            ++srcIndFracX.second;
            --srcIndFracX.first;
        }
        if (srcIndFracY.second < 0) {
            ++srcIndFracY.second;
            --srcIndFracY.first;
        }
 
        if (_srcGoodBBox.contains(lsst::geom::Point2I(srcIndFracX.first, srcIndFracY.first))) {
            // Offset source pixel index from kernel center to kernel corner (0, 0)
            // so we can convolveAtAPoint the pixels that overlap between source and kernel
            int srcStartX = srcIndFracX.first - _kernelCtr[0];
            int srcStartY = srcIndFracY.first - _kernelCtr[1];
 
            // Compute warped pixel
            double kSum = _setFracIndex(srcIndFracX.second, srcIndFracY.second);
 
            typename SrcImageT::const_xy_locator srcLoc = _srcImage.xy_at(srcStartX, srcStartY);
 
            *destXIter = lsst::afw::math::convolveAtAPoint<DestImageT, SrcImageT>(srcLoc, _xList, _yList);
            *destXIter *= relativeArea / kSum;
            return true;
        } else {
            // Edge pixel
            *destXIter = _padValue;
            return false;
        }
    }
 
    bool operator()(typename DestImageT::x_iterator &destXIter, lsst::geom::Point2D const &srcPos,
                    double relativeArea, lsst::afw::image::detail::MaskedImage_tag) {
        // Compute associated source pixel index as integer and nonnegative fractional parts;
        // the latter is used to compute the remapping kernel.
        std::pair<int, double> srcIndFracX = _srcImage.positionToIndex(srcPos[0], lsst::afw::image::X);
        std::pair<int, double> srcIndFracY = _srcImage.positionToIndex(srcPos[1], lsst::afw::image::Y);
        if (srcIndFracX.second < 0) {
            ++srcIndFracX.second;
            --srcIndFracX.first;
        }
        if (srcIndFracY.second < 0) {
            ++srcIndFracY.second;
            --srcIndFracY.first;
        }
 
        if (_srcGoodBBox.contains(lsst::geom::Point2I(srcIndFracX.first, srcIndFracY.first))) {
            // Offset source pixel index from kernel center to kernel corner (0, 0)
            // so we can convolveAtAPoint the pixels that overlap between source and kernel
            int srcStartX = srcIndFracX.first - _kernelCtr[0];
            int srcStartY = srcIndFracY.first - _kernelCtr[1];
 
            // Compute warped pixel
            double kSum = _setFracIndex(srcIndFracX.second, srcIndFracY.second);
 
            typename SrcImageT::const_xy_locator srcLoc = _srcImage.xy_at(srcStartX, srcStartY);
 
            *destXIter = lsst::afw::math::convolveAtAPoint<DestImageT, SrcImageT>(srcLoc, _xList, _yList);
            *destXIter *= relativeArea / kSum;
 
            if (_hasMaskKernel) {
                // compute mask value based on the mask kernel (replacing the value computed above)
                int maskStartX = srcIndFracX.first - _maskKernelCtr[0];
                int maskStartY = srcIndFracY.first - _maskKernelCtr[1];
 
                typename SrcImageT::Mask::const_xy_locator srcMaskLoc =
                        _srcImage.getMask()->xy_at(maskStartX, maskStartY);
 
                using k_iter = typename std::vector<lsst::afw::math::Kernel::Pixel>::const_iterator;
 
                typename DestImageT::Mask::SinglePixel destMaskValue = 0;
                for (double kValY : _maskYList) {
                    typename DestImageT::Mask::SinglePixel destMaskValueY = 0;
                    for (k_iter kernelXIter = _maskXList.begin(), xEnd = _maskXList.end();
                         kernelXIter != xEnd; ++kernelXIter, ++srcMaskLoc.x()) {
                        typename lsst::afw::math::Kernel::Pixel const kValX = *kernelXIter;
                        if (kValX != 0) {
                            destMaskValueY |= *srcMaskLoc;
                        }
                    }
 
                    if (kValY != 0) {
                        destMaskValue |= destMaskValueY;
                    }
 
                    srcMaskLoc += lsst::afw::image::detail::difference_type(-_maskXList.size(), 1);
                }
 
                destXIter.mask() = (destXIter.mask() & _growFullMask) | destMaskValue;
            }
            return true;
        } else {
            // Edge pixel
            *destXIter = _padValue;
            return false;
        }
    }
 
private:
    double _setFracIndex(double xFrac, double yFrac) {
        std::pair<double, double> srcFracInd(xFrac, yFrac);
        _kernelPtr->setKernelParameters(srcFracInd);
        double kSum = _kernelPtr->computeVectors(_xList, _yList, false);
        if (_maskKernelPtr) {
            _maskKernelPtr->setKernelParameters(srcFracInd);
            _maskKernelPtr->computeVectors(_maskXList, _maskYList, false);
        }
        return kSum;
    }
 
    SrcImageT _srcImage;
    std::shared_ptr<lsst::afw::math::SeparableKernel> _kernelPtr;
    std::shared_ptr<lsst::afw::math::SeparableKernel> _maskKernelPtr;
    bool _hasMaskKernel;
    lsst::geom::Point2I _kernelCtr;
    lsst::geom::Point2I _maskKernelCtr;
    lsst::afw::image::MaskPixel _growFullMask;
    std::vector<double> _xList;
    std::vector<double> _yList;
    std::vector<double> _maskXList;
    std::vector<double> _maskYList;
    typename DestImageT::SinglePixel _padValue;
    lsst::geom::Box2I const _srcGoodBBox;
};
}  // namespace detail
}  // namespace math
}  // namespace afw
}  // namespace lsst