lsst::afw::math::detail Namespace Reference

Classes
struct	ConvolveWithInterpolationWorkingImages
	kernel images used by convolveRegionWithInterpolation More...
class	KernelImagesForRegion
	A collection of Kernel images for special locations on a rectangular region of an image. More...
class	RowOfKernelImagesForRegion
	A row of KernelImagesForRegion. More...
class	SmoothedSpline
class	Spline
class	TautSpline
struct	TrapezoidalPacker
	A helper class ChebyshevBoundedField, for mapping trapezoidal matrices to 1-d arrays. More...
class	WarpAtOnePoint
	A functor that computes one warped pixel. More...

Functions
template<typename OutImageT , typename InImageT >
void	basicConvolve (OutImageT &convolvedImage, InImageT const &inImage, lsst::afw::math::Kernel const &kernel, lsst::afw::math::ConvolutionControl const &convolutionControl)
	Low-level convolution function that does not set edge pixels.
template<typename OutImageT , typename InImageT >
void	basicConvolve (OutImageT &convolvedImage, InImageT const &inImage, lsst::afw::math::DeltaFunctionKernel const &kernel, lsst::afw::math::ConvolutionControl const &convolutionControl)
	A version of basicConvolve that should be used when convolving delta function kernels.
template<typename OutImageT , typename InImageT >
void	basicConvolve (OutImageT &convolvedImage, InImageT const &inImage, lsst::afw::math::LinearCombinationKernel const &kernel, lsst::afw::math::ConvolutionControl const &convolutionControl)
	A version of basicConvolve that should be used when convolving a LinearCombinationKernel.
template<typename OutImageT , typename InImageT >
void	basicConvolve (OutImageT &convolvedImage, InImageT const &inImage, lsst::afw::math::SeparableKernel const &kernel, lsst::afw::math::ConvolutionControl const &convolutionControl)
	A version of basicConvolve that should be used when convolving separable kernels.
template<typename OutImageT , typename InImageT >
void	convolveWithBruteForce (OutImageT &convolvedImage, InImageT const &inImage, lsst::afw::math::Kernel const &kernel, lsst::afw::math::ConvolutionControl const &convolutionControl)
	Convolve an Image or MaskedImage with a Kernel by computing the kernel image at every point.
template<typename OutImageT , typename InImageT >
void	convolveWithInterpolation (OutImageT &outImage, InImageT const &inImage, lsst::afw::math::Kernel const &kernel, ConvolutionControl const &convolutionControl)
	Convolve an Image or MaskedImage with a spatially varying Kernel using linear interpolation.
template<typename OutImageT , typename InImageT >
void	convolveRegionWithInterpolation (OutImageT &outImage, InImageT const &inImage, KernelImagesForRegion const &region, ConvolveWithInterpolationWorkingImages &workingImages)
	Convolve a region of an Image or MaskedImage with a spatially varying Kernel using interpolation.
void	declareConvolve (lsst::utils::python::WrapperCollection &wrappers)
void	wrapConvolve (lsst::utils::python::WrapperCollection &wrappers)
void	wrapSpline (lsst::utils::python::WrapperCollection &)
	PYBIND11_MODULE (_detail, mod)

Detailed Description

Support code for lsst.afw.math

Function Documentation

basicConvolve() [1/4]

template<typename OutImageT , typename InImageT >

void lsst::afw::math::detail::basicConvolve	(	OutImageT &	convolvedImage,
		InImageT const &	inImage,
		lsst::afw::math::DeltaFunctionKernel const &	kernel,
		lsst::afw::math::ConvolutionControl const &	convolutionControl )

A version of basicConvolve that should be used when convolving delta function kernels.

Parameters

[out]	convolvedImage	convolved image
[in]	inImage	image to convolve
[in]	kernel	convolution kernel
[in]	convolutionControl	convolution control parameters

Definition at line 181 of file BasicConvolve.cc.

                                                                     {
    assert(!kernel.isSpatiallyVarying());
    assertDimensionsOK(convolvedImage, inImage, kernel);
 
    int const mImageWidth = inImage.getWidth();  // size of input region
    int const mImageHeight = inImage.getHeight();
    int const cnvWidth = mImageWidth + 1 - kernel.getWidth();
    int const cnvHeight = mImageHeight + 1 - kernel.getHeight();
    int const cnvStartX = kernel.getCtr().getX();
    int const cnvStartY = kernel.getCtr().getY();
    int const inStartX = kernel.getPixel().getX();
    int const inStartY = kernel.getPixel().getY();
 
    LOGL_DEBUG("TRACE2.lsst.afw.math.convolve.basicConvolve", "DeltaFunctionKernel basicConvolve");
 
    for (int i = 0; i < cnvHeight; ++i) {
        typename InImageT::x_iterator inPtr = inImage.x_at(inStartX, i + inStartY);
        for (typename OutImageT::x_iterator cnvPtr = convolvedImage.x_at(cnvStartX, i + cnvStartY),
                                            cnvEnd = cnvPtr + cnvWidth;
             cnvPtr != cnvEnd; ++cnvPtr, ++inPtr) {
            *cnvPtr = *inPtr;
        }
    }
}

basicConvolve() [2/4]

template<typename OutImageT , typename InImageT >

void lsst::afw::math::detail::basicConvolve	(	OutImageT &	convolvedImage,
		InImageT const &	inImage,
		lsst::afw::math::Kernel const &	kernel,
		lsst::afw::math::ConvolutionControl const &	convolutionControl )

Low-level convolution function that does not set edge pixels.

convolvedImage must be the same size as inImage. convolvedImage has a border in which the output pixels are not set. This border has size:

• kernel.getCtr().getX() along the left edge
• kernel.getCtr().getY() along the bottom edge
• kernel.getWidth() - 1 - kernel.getCtr().getX() along the right edge
• kernel.getHeight() - 1 - kernel.getCtr().getY() along the top edge

Parameters

[out]	convolvedImage	convolved image
[in]	inImage	image to convolve
[in]	kernel	convolution kernel
[in]	convolutionControl	convolution control parameters

Exceptions

lsst::pex::exceptions::InvalidParameterError	if convolvedImage dimensions != inImage dimensions
lsst::pex::exceptions::InvalidParameterError	if inImage smaller than kernel in width or height
lsst::pex::exceptions::InvalidParameterError	if kernel width or height < 1
std::bad_alloc	when allocation of CPU memory fails

Definition at line 141 of file BasicConvolve.cc.

                                                                     {
    // Because convolve isn't a method of Kernel we can't always use Kernel's vtbl to dynamically
    // dispatch the correct version of basicConvolve. The case that fails is convolving with a kernel
    // obtained from a pointer or reference to a Kernel (base class), e.g. as used in LinearCombinationKernel.
    if (IS_INSTANCE(kernel, math::DeltaFunctionKernel)) {
        LOGL_DEBUG("TRACE3.lsst.afw.math.convolve.basicConvolve",
                   "generic basicConvolve: dispatch to DeltaFunctionKernel basicConvolve");
        basicConvolve(convolvedImage, inImage, *dynamic_cast<math::DeltaFunctionKernel const*>(&kernel),
                      convolutionControl);
        return;
    } else if (IS_INSTANCE(kernel, math::SeparableKernel)) {
        LOGL_DEBUG("TRACE3.lsst.afw.math.convolve.basicConvolve",
                   "generic basicConvolve: dispatch to SeparableKernel basicConvolve");
        basicConvolve(convolvedImage, inImage, *dynamic_cast<math::SeparableKernel const*>(&kernel),
                      convolutionControl);
        return;
    } else if (IS_INSTANCE(kernel, math::LinearCombinationKernel) && kernel.isSpatiallyVarying()) {
        LOGL_DEBUG(
                "TRACE3.afw.math.convolve.basicConvolve",
                "generic basicConvolve: dispatch to spatially varying LinearCombinationKernel basicConvolve");
        basicConvolve(convolvedImage, inImage, *dynamic_cast<math::LinearCombinationKernel const*>(&kernel),
                      convolutionControl);
        return;
    }
    // OK, use general (and slower) form
    if (kernel.isSpatiallyVarying() && (convolutionControl.getMaxInterpolationDistance() > 1)) {
        // use linear interpolation
        LOGL_DEBUG("TRACE2.lsst.afw.math.convolve.basicConvolve",
                   "generic basicConvolve: using linear interpolation");
        convolveWithInterpolation(convolvedImage, inImage, kernel, convolutionControl);
 
    } else {
        // use brute force
        LOGL_DEBUG("TRACE2.lsst.afw.math.convolve.basicConvolve", "generic basicConvolve: using brute force");
        convolveWithBruteForce(convolvedImage, inImage, kernel, convolutionControl);
    }
}

basicConvolve() [3/4]

template<typename OutImageT , typename InImageT >

void lsst::afw::math::detail::basicConvolve	(	OutImageT &	convolvedImage,
		InImageT const &	inImage,
		lsst::afw::math::LinearCombinationKernel const &	kernel,
		lsst::afw::math::ConvolutionControl const &	convolutionControl )

A version of basicConvolve that should be used when convolving a LinearCombinationKernel.

The Algorithm:

• If the kernel is spatially varying and contains only DeltaFunctionKernels then convolves the input Image by each basis kernel in turn, solves the spatial model for that component and adds in the appropriate amount of the convolved image.
• In all other cases uses normal convolution

Parameters

[out]	convolvedImage	convolved image
[in]	inImage	image to convolve
[in]	kernel	convolution kernel
[in]	convolutionControl	convolution control parameters

Exceptions

lsst::pex::exceptions::InvalidParameterError	if convolvedImage dimensions != inImage dimensions
lsst::pex::exceptions::InvalidParameterError	if inImage smaller than kernel in width or height
lsst::pex::exceptions::InvalidParameterError	if kernel width or height < 1
std::bad_alloc	when allocation of CPU memory fails

Definition at line 209 of file BasicConvolve.cc.

                                                                     {
    if (!kernel.isSpatiallyVarying()) {
        // use the standard algorithm for the spatially invariant case
        LOGL_DEBUG("TRACE2.lsst.afw.math.convolve.basicConvolve",
                   "basicConvolve for LinearCombinationKernel: spatially invariant; using brute force");
        return convolveWithBruteForce(convolvedImage, inImage, kernel, convolutionControl.getDoNormalize());
    } else {
        // refactor the kernel if this is reasonable and possible;
        // then use the standard algorithm for the spatially varying case
        std::shared_ptr<Kernel> refKernelPtr;  // possibly refactored version of kernel
        if (static_cast<int>(kernel.getNKernelParameters()) > kernel.getNSpatialParameters()) {
            // refactoring will speed convolution, so try it
            refKernelPtr = kernel.refactor();
            if (!refKernelPtr) {
                refKernelPtr = kernel.clone();
            }
        } else {
            // too few basis kernels for refactoring to be worthwhile
            refKernelPtr = kernel.clone();
        }
        if (convolutionControl.getMaxInterpolationDistance() > 1) {
            LOGL_DEBUG("TRACE2.lsst.afw.math.convolve.basicConvolve",
                       "basicConvolve for LinearCombinationKernel: using interpolation");
            return convolveWithInterpolation(convolvedImage, inImage, *refKernelPtr, convolutionControl);
        } else {
            LOGL_DEBUG("TRACE2.lsst.afw.math.convolve.basicConvolve",
                       "basicConvolve for LinearCombinationKernel: maxInterpolationError < 0; using brute "
                       "force");
            return convolveWithBruteForce(convolvedImage, inImage, *refKernelPtr,
                                          convolutionControl.getDoNormalize());
        }
    }
}

basicConvolve() [4/4]

template<typename OutImageT , typename InImageT >

void lsst::afw::math::detail::basicConvolve	(	OutImageT &	convolvedImage,
		InImageT const &	inImage,
		lsst::afw::math::SeparableKernel const &	kernel,
		lsst::afw::math::ConvolutionControl const &	convolutionControl )

A version of basicConvolve that should be used when convolving separable kernels.

Parameters

[out]	convolvedImage	convolved image
[in]	inImage	image to convolve
[in]	kernel	convolution kernel
[in]	convolutionControl	convolution control parameters

Definition at line 246 of file BasicConvolve.cc.

                                                                     {
    using KernelPixel = typename math::Kernel::Pixel;
    using KernelVector = typename std::vector<KernelPixel>;
    using KernelIterator = KernelVector::const_iterator;
    using InXIterator = typename InImageT::const_x_iterator;
    using InXYLocator = typename InImageT::const_xy_locator;
    using OutXIterator = typename OutImageT::x_iterator;
    using OutYIterator = typename OutImageT::y_iterator;
    using OutPixel = typename OutImageT::SinglePixel;
 
    assertDimensionsOK(convolvedImage, inImage, kernel);
 
    lsst::geom::Box2I const fullBBox = inImage.getBBox(image::LOCAL);
    lsst::geom::Box2I const goodBBox = kernel.shrinkBBox(fullBBox);
 
    KernelVector kernelXVec(kernel.getWidth());
    KernelVector kernelYVec(kernel.getHeight());
 
    if (kernel.isSpatiallyVarying()) {
        LOGL_DEBUG("TRACE2.lsst.afw.math.convolve.basicConvolve",
                   "SeparableKernel basicConvolve: kernel is spatially varying");
 
        for (int cnvY = goodBBox.getMinY(); cnvY <= goodBBox.getMaxY(); ++cnvY) {
            double const rowPos = inImage.indexToPosition(cnvY, image::Y);
 
            InXYLocator inImLoc = inImage.xy_at(0, cnvY - goodBBox.getMinY());
            OutXIterator cnvXIter = convolvedImage.row_begin(cnvY) + goodBBox.getMinX();
            for (int cnvX = goodBBox.getMinX(); cnvX <= goodBBox.getMaxX();
                 ++cnvX, ++inImLoc.x(), ++cnvXIter) {
                double const colPos = inImage.indexToPosition(cnvX, image::X);
 
                KernelPixel kSum = kernel.computeVectors(kernelXVec, kernelYVec,
                                                         convolutionControl.getDoNormalize(), colPos, rowPos);
 
                // why does this trigger warnings? It did not in the past.
                *cnvXIter = math::convolveAtAPoint<OutImageT, InImageT>(inImLoc, kernelXVec, kernelYVec);
                if (convolutionControl.getDoNormalize()) {
                    *cnvXIter = *cnvXIter / kSum;
                }
            }
        }
    } else {
        // kernel is spatially invariant
        // The basic sequence:
        // - For each output row:
        // - Compute x-convolved data: a kernel height's strip of input image convolved with kernel x vector
        // - Compute one row of output by dotting each column of x-convolved data with the kernel y vector
        // The x-convolved data is stored in a kernel-height by good-width buffer.
        // This is circular buffer along y (to avoid shifting pixels before setting each new row);
        // so for each new row the kernel y vector is rotated to match the order of the x-convolved data.
 
        LOGL_DEBUG("TRACE2.lsst.afw.math.convolve.basicConvolve",
                   "SeparableKernel basicConvolve: kernel is spatially invariant");
 
        kernel.computeVectors(kernelXVec, kernelYVec, convolutionControl.getDoNormalize());
        KernelIterator const kernelXVecBegin = kernelXVec.begin();
        KernelIterator const kernelYVecBegin = kernelYVec.begin();
 
        // buffer for x-convolved data
        OutImageT buffer(lsst::geom::Extent2I(goodBBox.getWidth(), kernel.getHeight()));
 
        // pre-fill x-convolved data buffer with all but one row of data
        int yInd = 0;  // during initial fill bufY = inImageY
        int const yPrefillEnd = buffer.getHeight() - 1;
        for (; yInd < yPrefillEnd; ++yInd) {
            OutXIterator bufXIter = buffer.x_at(0, yInd);
            OutXIterator const bufXEnd = buffer.x_at(goodBBox.getWidth(), yInd);
            InXIterator inXIter = inImage.x_at(0, yInd);
            for (; bufXIter != bufXEnd; ++bufXIter, ++inXIter) {
                *bufXIter = kernelDotProduct<OutPixel, InXIterator, KernelIterator, KernelPixel>(
                        inXIter, kernelXVecBegin, kernel.getWidth());
            }
        }
 
        // compute output pixels using the sequence described above
        int inY = yPrefillEnd;
        int bufY = yPrefillEnd;
        int cnvY = goodBBox.getMinY();
        while (true) {
            // fill next buffer row and compute output row
            InXIterator inXIter = inImage.x_at(0, inY);
            OutXIterator bufXIter = buffer.x_at(0, bufY);
            OutXIterator cnvXIter = convolvedImage.x_at(goodBBox.getMinX(), cnvY);
            for (int bufX = 0; bufX < goodBBox.getWidth(); ++bufX, ++cnvXIter, ++bufXIter, ++inXIter) {
                // note: bufXIter points to the row of the buffer that is being updated,
                // whereas bufYIter points to row 0 of the buffer
                *bufXIter = kernelDotProduct<OutPixel, InXIterator, KernelIterator, KernelPixel>(
                        inXIter, kernelXVecBegin, kernel.getWidth());
 
                OutYIterator bufYIter = buffer.y_at(bufX, 0);
                *cnvXIter = kernelDotProduct<OutPixel, OutYIterator, KernelIterator, KernelPixel>(
                        bufYIter, kernelYVecBegin, kernel.getHeight());
            }
 
            // test for done now, instead of the start of the loop,
            // to avoid an unnecessary extra rotation of the kernel Y vector
            if (cnvY >= goodBBox.getMaxY()) break;
 
            // update y indices, including bufY, and rotate the kernel y vector to match
            ++inY;
            bufY = (bufY + 1) % kernel.getHeight();
            ++cnvY;
            std::rotate(kernelYVec.begin(), kernelYVec.end() - 1, kernelYVec.end());
        }
    }
}

convolveRegionWithInterpolation()

template<typename OutImageT , typename InImageT >

void lsst::afw::math::detail::convolveRegionWithInterpolation	(	OutImageT &	outImage,
		InImageT const &	inImage,
		KernelImagesForRegion const &	region,
		ConvolveWithInterpolationWorkingImages &	workingImages )

Convolve a region of an Image or MaskedImage with a spatially varying Kernel using interpolation.

This is a low-level convolution function that does not set edge pixels.

Parameters

[out]	outImage	convolved image = inImage convolved with kernel
[in]	inImage	input image
[in]	region	kernel image region over which to convolve
[in]	workingImages	working kernel images

Warning

: this is a low-level routine that performs no bounds checking.

Definition at line 88 of file ConvolveWithInterpolation.cc.

                                                                                            {
    using OutLocator = typename OutImageT::xy_locator;
    using InConstLocator = typename InImageT::const_xy_locator;
    using KernelImage = KernelImagesForRegion::Image;
    using KernelConstLocator = KernelImage::const_xy_locator;
 
    std::shared_ptr<Kernel const> kernelPtr = region.getKernel();
    lsst::geom::Extent2I const kernelDimensions(kernelPtr->getDimensions());
    workingImages.leftImage.assign(*region.getImage(KernelImagesForRegion::BOTTOM_LEFT));
    workingImages.rightImage.assign(*region.getImage(KernelImagesForRegion::BOTTOM_RIGHT));
    workingImages.kernelImage.assign(workingImages.leftImage);
 
    lsst::geom::Box2I const goodBBox = region.getBBox();
    lsst::geom::Box2I const fullBBox = kernelPtr->growBBox(goodBBox);
 
    // top and right images are computed one beyond bbox boundary,
    // so the distance between edge images is bbox width/height pixels
    double xfrac = 1.0 / static_cast<double>(goodBBox.getWidth());
    double yfrac = 1.0 / static_cast<double>(goodBBox.getHeight());
    math::scaledPlus(workingImages.leftDeltaImage, yfrac, *region.getImage(KernelImagesForRegion::TOP_LEFT),
                     -yfrac, workingImages.leftImage);
    math::scaledPlus(workingImages.rightDeltaImage, yfrac, *region.getImage(KernelImagesForRegion::TOP_RIGHT),
                     -yfrac, workingImages.rightImage);
 
    KernelConstLocator const kernelLocator = workingImages.kernelImage.xy_at(0, 0);
 
    // The loop is a bit odd for efficiency: the initial value of workingImages.kernelImage
    // and related kernel images are set when they are allocated,
    // so they are not computed in the loop until after the convolution; to save cpu cycles
    // they are not computed at all for the last iteration.
    InConstLocator inLocator = inImage.xy_at(fullBBox.getMinX(), fullBBox.getMinY());
    OutLocator outLocator = outImage.xy_at(goodBBox.getMinX(), goodBBox.getMinY());
    for (int j = 0;;) {
        auto inLocatorInitialPosition = inLocator;
        auto outLocatorInitialPosition = outLocator;
        math::scaledPlus(workingImages.deltaImage, xfrac, workingImages.rightImage, -xfrac,
                         workingImages.leftImage);
        for (int i = 0;;) {
            *outLocator = math::convolveAtAPoint<OutImageT, InImageT>(
                    inLocator, kernelLocator, kernelDimensions.getX(), kernelDimensions.getY());
            ++outLocator.x();
            ++inLocator.x();
            ++i;
            if (i >= goodBBox.getWidth()) {
                break;
            }
            workingImages.kernelImage += workingImages.deltaImage;
        }
 
        ++j;
        if (j >= goodBBox.getHeight()) {
            break;
        }
        workingImages.leftImage += workingImages.leftDeltaImage;
        workingImages.rightImage += workingImages.rightDeltaImage;
        workingImages.kernelImage.assign(workingImages.leftImage);
 
        // Workaround for DM-5822
        // Boost GIL locator won't decrement in x-dimension for some strange and still
        // not understood reason. So instead store position at start of previous row,
        // reset and move down.
        inLocator = inLocatorInitialPosition;
        outLocator = outLocatorInitialPosition;
        ++inLocator.y();
        ++outLocator.y();
    }
}

convolveWithBruteForce()

template<typename OutImageT , typename InImageT >

void lsst::afw::math::detail::convolveWithBruteForce	(	OutImageT &	convolvedImage,
		InImageT const &	inImage,
		lsst::afw::math::Kernel const &	kernel,
		lsst::afw::math::ConvolutionControl const &	convolutionControl )

Convolve an Image or MaskedImage with a Kernel by computing the kernel image at every point.

(If the kernel is not spatially varying then only compute it once).

convolvedImage must be the same size as inImage. convolvedImage has a border in which the output pixels are not set. This border has size:

• kernel.getCtr().getX() along the left edge
• kernel.getCtr().getY() along the bottom edge
• kernel.getWidth() - 1 - kernel.getCtr().getX() along the right edge
• kernel.getHeight() - 1 - kernel.getCtr().getY() along the top edge

Parameters

[out]	convolvedImage	convolved image
[in]	inImage	image to convolve
[in]	kernel	convolution kernel
[in]	convolutionControl	convolution control parameters

Exceptions

lsst::pex::exceptions::InvalidParameterError	if convolvedImage dimensions != inImage dimensions
lsst::pex::exceptions::InvalidParameterError	if inImage smaller than kernel in width or height
lsst::pex::exceptions::InvalidParameterError	if kernel width or height < 1
std::bad_alloc	when allocation of CPU memory fails

Warning

Low-level convolution function that does not set edge pixels.

Definition at line 355 of file BasicConvolve.cc.

                                                                              {
    bool doNormalize = convolutionControl.getDoNormalize();
 
    using KernelPixel = typename math::Kernel::Pixel;
    using KernelImage = image::Image<KernelPixel>;
 
    using KernelXIterator = typename KernelImage::const_x_iterator;
    using KernelXYLocator = typename KernelImage::const_xy_locator;
    using InXIterator = typename InImageT::const_x_iterator;
    using InXYLocator = typename InImageT::const_xy_locator;
    using OutXIterator = typename OutImageT::x_iterator;
    using OutPixel = typename OutImageT::SinglePixel;
 
    assertDimensionsOK(convolvedImage, inImage, kernel);
 
    int const inImageWidth = inImage.getWidth();
    int const inImageHeight = inImage.getHeight();
    int const kWidth = kernel.getWidth();
    int const kHeight = kernel.getHeight();
    int const cnvWidth = inImageWidth + 1 - kernel.getWidth();
    int const cnvHeight = inImageHeight + 1 - kernel.getHeight();
    int const cnvStartX = kernel.getCtr().getX();
    int const cnvStartY = kernel.getCtr().getY();
    int const cnvEndX = cnvStartX + cnvWidth;   // end index + 1
    int const cnvEndY = cnvStartY + cnvHeight;  // end index + 1
 
    KernelImage kernelImage(kernel.getDimensions());
    KernelXYLocator const kernelLoc = kernelImage.xy_at(0, 0);
 
    if (kernel.isSpatiallyVarying()) {
        LOGL_DEBUG("TRACE4.lsst.afw.math.convolve.convolveWithBruteForce",
                   "convolveWithBruteForce: kernel is spatially varying");
 
        for (int cnvY = cnvStartY; cnvY != cnvEndY; ++cnvY) {
            double const rowPos = inImage.indexToPosition(cnvY, image::Y);
 
            InXYLocator inImLoc = inImage.xy_at(0, cnvY - cnvStartY);
            OutXIterator cnvXIter = convolvedImage.x_at(cnvStartX, cnvY);
            for (int cnvX = cnvStartX; cnvX != cnvEndX; ++cnvX, ++inImLoc.x(), ++cnvXIter) {
                double const colPos = inImage.indexToPosition(cnvX, image::X);
 
                KernelPixel kSum = kernel.computeImage(kernelImage, false, colPos, rowPos);
                *cnvXIter = math::convolveAtAPoint<OutImageT, InImageT>(inImLoc, kernelLoc, kWidth, kHeight);
                if (doNormalize) {
                    *cnvXIter = *cnvXIter / kSum;
                }
            }
        }
    } else {
        LOGL_DEBUG("TRACE4.lsst.afw.math.convolve.convolveWithBruteForce",
                   "convolveWithBruteForce: kernel is spatially invariant");
 
        (void)kernel.computeImage(kernelImage, doNormalize);
 
        for (int inStartY = 0, cnvY = cnvStartY; inStartY < cnvHeight; ++inStartY, ++cnvY) {
            KernelXIterator kernelXIter = kernelImage.x_at(0, 0);
            InXIterator inXIter = inImage.x_at(0, inStartY);
            OutXIterator cnvXIter = convolvedImage.x_at(cnvStartX, cnvY);
            for (int x = 0; x < cnvWidth; ++x, ++cnvXIter, ++inXIter) {
                *cnvXIter = kernelDotProduct<OutPixel, InXIterator, KernelXIterator, KernelPixel>(
                        inXIter, kernelXIter, kWidth);
            }
            for (int kernelY = 1, inY = inStartY + 1; kernelY < kHeight; ++inY, ++kernelY) {
                KernelXIterator kernelXIter = kernelImage.x_at(0, kernelY);
                InXIterator inXIter = inImage.x_at(0, inY);
                OutXIterator cnvXIter = convolvedImage.x_at(cnvStartX, cnvY);
                for (int x = 0; x < cnvWidth; ++x, ++cnvXIter, ++inXIter) {
                    *cnvXIter += kernelDotProduct<OutPixel, InXIterator, KernelXIterator, KernelPixel>(
                            inXIter, kernelXIter, kWidth);
                }
            }
        }
    }
}

convolveWithInterpolation()

template<typename OutImageT , typename InImageT >

void lsst::afw::math::detail::convolveWithInterpolation	(	OutImageT &	outImage,
		InImageT const &	inImage,
		lsst::afw::math::Kernel const &	kernel,
		math::ConvolutionControl const &	convolutionControl )

Convolve an Image or MaskedImage with a spatially varying Kernel using linear interpolation.

This is a low-level convolution function that does not set edge pixels.

The algorithm is as follows:

• divide the image into regions whose size is no larger than maxInterpolationDistance
• for each region:
  • convolve it using convolveRegionWithInterpolation (which see)

Note that this routine will also work with spatially invariant kernels, but not efficiently.

Parameters

[out]	outImage	convolved image = inImage convolved with kernel
[in]	inImage	input image
[in]	kernel	convolution kernel
[in]	convolutionControl	convolution control parameters

Exceptions

lsst::pex::exceptions::InvalidParameterError

if outImage is not the same size as inImage

Definition at line 45 of file ConvolveWithInterpolation.cc.

                                                                                 {
    if (outImage.getDimensions() != inImage.getDimensions()) {
        std::ostringstream os;
        os << "outImage dimensions = ( " << outImage.getWidth() << ", " << outImage.getHeight() << ") != ("
           << inImage.getWidth() << ", " << inImage.getHeight() << ") = inImage dimensions";
        throw LSST_EXCEPT(pexExcept::InvalidParameterError, os.str());
    }
 
    // compute region covering good area of output image
    lsst::geom::Box2I fullBBox = lsst::geom::Box2I(
            lsst::geom::Point2I(0, 0), lsst::geom::Extent2I(outImage.getWidth(), outImage.getHeight()));
    lsst::geom::Box2I goodBBox = kernel.shrinkBBox(fullBBox);
    KernelImagesForRegion goodRegion(KernelImagesForRegion(kernel.clone(), goodBBox, inImage.getXY0(),
                                                           convolutionControl.getDoNormalize()));
    LOGL_DEBUG("TRACE5.lsst.afw.math.convolve.convolveWithInterpolation",
               "convolveWithInterpolation: full bbox minimum=(%d, %d), extent=(%d, %d)", fullBBox.getMinX(),
               fullBBox.getMinY(), fullBBox.getWidth(), fullBBox.getHeight());
    LOGL_DEBUG("TRACE5.lsst.afw.math.convolve.convolveWithInterpolation",
               "convolveWithInterpolation: goodRegion bbox minimum=(%d, %d), extent=(%d, %d)",
               goodRegion.getBBox().getMinX(), goodRegion.getBBox().getMinY(),
               goodRegion.getBBox().getWidth(), goodRegion.getBBox().getHeight());
 
    // divide good region into subregions small enough to interpolate over
    int nx = 1 + (goodBBox.getWidth() / convolutionControl.getMaxInterpolationDistance());
    int ny = 1 + (goodBBox.getHeight() / convolutionControl.getMaxInterpolationDistance());
    LOGL_DEBUG("TRACE3.lsst.afw.math.convolve.convolveWithInterpolation",
               "convolveWithInterpolation: divide into %d x %d subregions", nx, ny);
 
    ConvolveWithInterpolationWorkingImages workingImages(kernel.getDimensions());
    RowOfKernelImagesForRegion regionRow(nx, ny);
    while (goodRegion.computeNextRow(regionRow)) {
        for (auto const &rgnIter : regionRow) {
            LOGL_DEBUG("TRACE5.lsst.afw.math.convolve.convolveWithInterpolation",
                       "convolveWithInterpolation: bbox minimum=(%d, %d), extent=(%d, %d)",
                       rgnIter->getBBox().getMinX(), rgnIter->getBBox().getMinY(),
                       rgnIter->getBBox().getWidth(), rgnIter->getBBox().getHeight());
            convolveRegionWithInterpolation(outImage, inImage, *rgnIter, workingImages);
        }
    }
}

declareConvolve()

void lsst::afw::math::detail::declareConvolve

(

lsst::utils::python::WrapperCollection &

wrappers

)

Definition at line 69 of file _convolve.cc.

                                                                 {
    using PyClass = py::class_<KernelImagesForRegion, std::shared_ptr<KernelImagesForRegion>>;
    auto clsKernelImagesForRegion =
            wrappers.wrapType(PyClass(wrappers.module, "KernelImagesForRegion"), [](auto &mod, auto &cls) {
                cls.def(py::init<KernelImagesForRegion::KernelConstPtr, lsst::geom::Box2I const &,
                                 lsst::geom::Point2I const &, bool>(),
                        "kernelPtr"_a, "bbox"_a, "xy0"_a, "doNormalize"_a);
                cls.def(py::init<KernelImagesForRegion::KernelConstPtr, lsst::geom::Box2I const &,
                                 lsst::geom::Point2I const &, bool, KernelImagesForRegion::ImagePtr,
                                 KernelImagesForRegion::ImagePtr, KernelImagesForRegion::ImagePtr,
                                 KernelImagesForRegion::ImagePtr>(),
                        "kernelPtr"_a, "bbox"_a, "xy0"_a, "doNormalize"_a, "bottomLeftImagePtr"_a,
                        "bottomRightImagePtr"_a, "topLeftImagePtr"_a, "topRightImagePtr"_a);
 
                cls.def("getBBox", &KernelImagesForRegion::getBBox);
                cls.def("getXY0", &KernelImagesForRegion::getXY0);
                cls.def("getDoNormalize", &KernelImagesForRegion::getDoNormalize);
                cls.def("getImage", &KernelImagesForRegion::getImage);
                cls.def("getKernel", &KernelImagesForRegion::getKernel);
                cls.def("getPixelIndex", &KernelImagesForRegion::getPixelIndex);
                cls.def("computeNextRow", &KernelImagesForRegion::computeNextRow);
                cls.def_static("getMinInterpolationSize", KernelImagesForRegion::getMinInterpolationSize);
            });
 
    wrappers.wrapType(py::enum_<KernelImagesForRegion::Location>(clsKernelImagesForRegion, "Location"),
                      [](auto &mod, auto &enm) {
                          enm.value("BOTTOM_LEFT", KernelImagesForRegion::Location::BOTTOM_LEFT);
                          enm.value("BOTTOM_RIGHT", KernelImagesForRegion::Location::BOTTOM_RIGHT);
                          enm.value("TOP_LEFT", KernelImagesForRegion::Location::TOP_LEFT);
                          enm.value("TOP_RIGHT", KernelImagesForRegion::Location::TOP_RIGHT);
                          enm.export_values();
                      });
 
    wrappers.wrapType(py::class_<RowOfKernelImagesForRegion, std::shared_ptr<RowOfKernelImagesForRegion>>(
                              wrappers.module, "RowOfKernelImagesForRegion"),
                      [](auto &mod, auto &cls) {
                          cls.def(py::init<int, int>(), "nx"_a, "ny"_a);
 
                          cls.def("front", &RowOfKernelImagesForRegion::front);
                          cls.def("back", &RowOfKernelImagesForRegion::back);
                          cls.def("getNX", &RowOfKernelImagesForRegion::getNX);
                          cls.def("getNY", &RowOfKernelImagesForRegion::getNY);
                          cls.def("getYInd", &RowOfKernelImagesForRegion::getYInd);
                          cls.def("getRegion", &RowOfKernelImagesForRegion::getRegion);
                          cls.def("hasData", &RowOfKernelImagesForRegion::hasData);
                          cls.def("isLastRow", &RowOfKernelImagesForRegion::isLastRow);
                          cls.def("incrYInd", &RowOfKernelImagesForRegion::incrYInd);
                      });
}

PYBIND11_MODULE()

lsst::afw::math::detail::PYBIND11_MODULE	(	_detail	,
		mod	)

Definition at line 11 of file _detail.cc.

                              {
    lsst::utils::python::WrapperCollection wrappers(mod, "lsst.afw.math.detail");
    wrapConvolve(wrappers);
    wrapSpline(wrappers);
    wrappers.finish();
}

wrapConvolve()

void lsst::afw::math::detail::wrapConvolve

(

lsst::utils::python::WrapperCollection &

wrappers

)

Definition at line 118 of file _convolve.cc.

                                                              {
    wrappers.addSignatureDependency("lsst.afw.image");
    wrappers.addSignatureDependency("lsst.afw.math");
    declareConvolve(wrappers);
    declareAll<double, double>(wrappers);
    declareAll<double, float>(wrappers);
    declareAll<double, int>(wrappers);
    declareAll<double, std::uint16_t>(wrappers);
    declareAll<float, float>(wrappers);
    declareAll<float, int>(wrappers);
    declareAll<float, std::uint16_t>(wrappers);
    declareAll<int, int>(wrappers);
    declareAll<std::uint16_t, std::uint16_t>(wrappers);
}

wrapSpline()

void lsst::afw::math::detail::wrapSpline

(

lsst::utils::python::WrapperCollection &

wrappers

)

Definition at line 38 of file _spline.cc.

                                                            {
    /* Module level */
    wrappers.wrapType(py::class_<Spline>(wrappers.module, "Spline"), [](auto &mod, auto &cls) {
        cls.def("interpolate", &Spline::interpolate);
        cls.def("derivative", &Spline::derivative);
    });
    auto clsTautSpline = wrappers.wrapType(
            py::class_<TautSpline, Spline>(wrappers.module, "TautSpline"), [](auto &mod, auto &cls) {
                cls.def(py::init<std::vector<double> const &, std::vector<double> const &, double const,
                                 TautSpline::Symmetry>(),
                        "x"_a, "y"_a, "gamma"_a = 0, "type"_a = lsst::afw::math::detail::TautSpline::Unknown);
                cls.def("roots", &TautSpline::roots);
            });
    /* Member types and enums */
    wrappers.wrapType(py::enum_<TautSpline::Symmetry>(clsTautSpline, "Symmetry"), [](auto &mod, auto &enm) {
        enm.value("Unknown", TautSpline::Symmetry::Unknown);
        enm.value("Odd", TautSpline::Symmetry::Odd);
        enm.value("Even", TautSpline::Symmetry::Even);
        enm.export_values();
    });
}