LSSTApplications  18.1.0
LSSTDataManagementBasePackage
ConvolveWithInterpolation.cc
Go to the documentation of this file.
1 // -*- LSST-C++ -*-
2 
3 /*
4  * LSST Data Management System
5  * Copyright 2008, 2009, 2010 LSST Corporation.
6  *
7  * This product includes software developed by the
8  * LSST Project (http://www.lsst.org/).
9  *
10  * This program is free software: you can redistribute it and/or modify
11  * it under the terms of the GNU General Public License as published by
12  * the Free Software Foundation, either version 3 of the License, or
13  * (at your option) any later version.
14  *
15  * This program is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  * GNU General Public License for more details.
19  *
20  * You should have received a copy of the LSST License Statement and
21  * the GNU General Public License along with this program. If not,
22  * see <http://www.lsstcorp.org/LegalNotices/>.
23  */
24 
25 /*
26  * Definition of convolveWithInterpolation and helper functions declared in detail/ConvolveImage.h
27  */
28 #include <algorithm>
29 #include <cmath>
30 #include <cstdint>
31 #include <sstream>
32 #include <vector>
33 #include <iostream>
34 
35 #include "lsst/pex/exceptions.h"
36 #include "lsst/log/Log.h"
37 #include "lsst/geom.h"
39 #include "lsst/afw/math/Kernel.h"
41 
43 
44 namespace lsst {
45 namespace afw {
46 namespace math {
47 namespace detail {
48 
49 template <typename OutImageT, typename InImageT>
50 void convolveWithInterpolation(OutImageT &outImage, InImageT const &inImage, math::Kernel const &kernel,
51  math::ConvolutionControl const &convolutionControl) {
52  if (outImage.getDimensions() != inImage.getDimensions()) {
54  os << "outImage dimensions = ( " << outImage.getWidth() << ", " << outImage.getHeight() << ") != ("
55  << inImage.getWidth() << ", " << inImage.getHeight() << ") = inImage dimensions";
57  }
58 
59  // compute region covering good area of output image
61  lsst::geom::Point2I(0, 0), lsst::geom::Extent2I(outImage.getWidth(), outImage.getHeight()));
62  lsst::geom::Box2I goodBBox = kernel.shrinkBBox(fullBBox);
63  KernelImagesForRegion goodRegion(KernelImagesForRegion(kernel.clone(), goodBBox, inImage.getXY0(),
64  convolutionControl.getDoNormalize()));
65  LOGL_DEBUG("TRACE5.afw.math.convolve.convolveWithInterpolation",
66  "convolveWithInterpolation: full bbox minimum=(%d, %d), extent=(%d, %d)", fullBBox.getMinX(),
67  fullBBox.getMinY(), fullBBox.getWidth(), fullBBox.getHeight());
68  LOGL_DEBUG("TRACE5.afw.math.convolve.convolveWithInterpolation",
69  "convolveWithInterpolation: goodRegion bbox minimum=(%d, %d), extent=(%d, %d)",
70  goodRegion.getBBox().getMinX(), goodRegion.getBBox().getMinY(),
71  goodRegion.getBBox().getWidth(), goodRegion.getBBox().getHeight());
72 
73  // divide good region into subregions small enough to interpolate over
74  int nx = 1 + (goodBBox.getWidth() / convolutionControl.getMaxInterpolationDistance());
75  int ny = 1 + (goodBBox.getHeight() / convolutionControl.getMaxInterpolationDistance());
76  LOGL_DEBUG("TRACE3.afw.math.convolve.convolveWithInterpolation",
77  "convolveWithInterpolation: divide into %d x %d subregions", nx, ny);
78 
80  RowOfKernelImagesForRegion regionRow(nx, ny);
81  while (goodRegion.computeNextRow(regionRow)) {
82  for (RowOfKernelImagesForRegion::ConstIterator rgnIter = regionRow.begin(), rgnEnd = regionRow.end();
83  rgnIter != rgnEnd; ++rgnIter) {
84  LOGL_DEBUG("TRACE5.afw.math.convolve.convolveWithInterpolation",
85  "convolveWithInterpolation: bbox minimum=(%d, %d), extent=(%d, %d)",
86  (*rgnIter)->getBBox().getMinX(), (*rgnIter)->getBBox().getMinY(),
87  (*rgnIter)->getBBox().getWidth(), (*rgnIter)->getBBox().getHeight());
88  convolveRegionWithInterpolation(outImage, inImage, **rgnIter, workingImages);
89  }
90  }
91 }
92 
93 template <typename OutImageT, typename InImageT>
94 void convolveRegionWithInterpolation(OutImageT &outImage, InImageT const &inImage,
97  typedef typename OutImageT::xy_locator OutLocator;
98  typedef typename InImageT::const_xy_locator InConstLocator;
99  typedef KernelImagesForRegion::Image KernelImage;
100  typedef KernelImage::const_xy_locator KernelConstLocator;
101 
102  std::shared_ptr<Kernel const> kernelPtr = region.getKernel();
103  lsst::geom::Extent2I const kernelDimensions(kernelPtr->getDimensions());
106  workingImages.kernelImage.assign(workingImages.leftImage);
107 
108  lsst::geom::Box2I const goodBBox = region.getBBox();
109  lsst::geom::Box2I const fullBBox = kernelPtr->growBBox(goodBBox);
110 
111  // top and right images are computed one beyond bbox boundary,
112  // so the distance between edge images is bbox width/height pixels
113  double xfrac = 1.0 / static_cast<double>(goodBBox.getWidth());
114  double yfrac = 1.0 / static_cast<double>(goodBBox.getHeight());
116  -yfrac, workingImages.leftImage);
118  -yfrac, workingImages.rightImage);
119 
120  KernelConstLocator const kernelLocator = workingImages.kernelImage.xy_at(0, 0);
121 
122  // The loop is a bit odd for efficiency: the initial value of workingImages.kernelImage
123  // and related kernel images are set when they are allocated,
124  // so they are not computed in the loop until after the convolution; to save cpu cycles
125  // they are not computed at all for the last iteration.
126  InConstLocator inLocator = inImage.xy_at(fullBBox.getMinX(), fullBBox.getMinY());
127  OutLocator outLocator = outImage.xy_at(goodBBox.getMinX(), goodBBox.getMinY());
128  for (int j = 0;;) {
129  auto inLocatorInitialPosition = inLocator;
130  auto outLocatorInitialPosition = outLocator;
131  math::scaledPlus(workingImages.deltaImage, xfrac, workingImages.rightImage, -xfrac,
132  workingImages.leftImage);
133  for (int i = 0;;) {
134  *outLocator = math::convolveAtAPoint<OutImageT, InImageT>(
135  inLocator, kernelLocator, kernelDimensions.getX(), kernelDimensions.getY());
136  ++outLocator.x();
137  ++inLocator.x();
138  ++i;
139  if (i >= goodBBox.getWidth()) {
140  break;
141  }
142  workingImages.kernelImage += workingImages.deltaImage;
143  }
144 
145  ++j;
146  if (j >= goodBBox.getHeight()) {
147  break;
148  }
149  workingImages.leftImage += workingImages.leftDeltaImage;
150  workingImages.rightImage += workingImages.rightDeltaImage;
151  workingImages.kernelImage.assign(workingImages.leftImage);
152 
153  // Workaround for DM-5822
154  // Boost GIL locator won't decrement in x-dimension for some strange and still
155  // not understood reason. So instead store position at start of previous row,
156  // reset and move down.
157  inLocator = inLocatorInitialPosition;
158  outLocator = outLocatorInitialPosition;
159  ++inLocator.y();
160  ++outLocator.y();
161  }
162 }
163 
164 /*
165  * Explicit instantiation
166  */
168 #define IMAGE(PIXTYPE) image::Image<PIXTYPE>
169 #define MASKEDIMAGE(PIXTYPE) image::MaskedImage<PIXTYPE, image::MaskPixel, image::VariancePixel>
170 #define NL /* */
171 // Instantiate Image or MaskedImage versions
172 #define INSTANTIATE_IM_OR_MI(IMGMACRO, OUTPIXTYPE, INPIXTYPE) \
173  template void convolveWithInterpolation(IMGMACRO(OUTPIXTYPE) &, IMGMACRO(INPIXTYPE) const &, \
174  math::Kernel const &, math::ConvolutionControl const &); \
175  NL template void convolveRegionWithInterpolation(IMGMACRO(OUTPIXTYPE) &, IMGMACRO(INPIXTYPE) const &, \
176  KernelImagesForRegion const &, \
177  ConvolveWithInterpolationWorkingImages &);
178 // Instantiate both Image and MaskedImage versions
179 #define INSTANTIATE(OUTPIXTYPE, INPIXTYPE) \
180  INSTANTIATE_IM_OR_MI(IMAGE, OUTPIXTYPE, INPIXTYPE) \
181  INSTANTIATE_IM_OR_MI(MASKEDIMAGE, OUTPIXTYPE, INPIXTYPE)
182 
183 INSTANTIATE(double, double)
184 INSTANTIATE(double, float)
185 INSTANTIATE(double, int)
186 INSTANTIATE(double, std::uint16_t)
187 INSTANTIATE(float, float)
188 INSTANTIATE(float, int)
190 INSTANTIATE(int, int)
193 } // namespace detail
194 } // namespace math
195 } // namespace afw
196 } // namespace lsst
void scaledPlus(OutImageT &outImage, double c1, InImageT const &inImage1, double c2, InImageT const &inImage2)
Compute the scaled sum of two images.
lsst::geom::Box2I getBBox() const
Get the bounding box for the region.
Definition: Convolve.h:234
int getHeight() const noexcept
Definition: Box.h:175
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...
Parameters to control convolution.
Definition: ConvolveImage.h:50
ImagePtr getImage(Location location) const
Return the image and sum at the specified location.
#define LOGL_DEBUG(logger, message...)
Log a debug-level message using a varargs/printf style interface.
Definition: Log.h:489
kernel images used by convolveRegionWithInterpolation
Definition: Convolve.h:418
LSST DM logging module built on log4cxx.
A base class for image defects.
bool computeNextRow(RowOfKernelImagesForRegion &regionRow) const
Compute next row of subregions.
void assign(ImageBase const &rhs, lsst::geom::Box2I const &bbox=lsst::geom::Box2I(), ImageOrigin origin=PARENT)
Copy pixels from another image to a specified subregion of this image.
Definition: Image.cc:177
T str(T... args)
int getWidth() const noexcept
Definition: Box.h:174
A row of KernelImagesForRegion.
Definition: Convolve.h:334
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.
lsst::geom::Extent2I const getDimensions() const
Return the Kernel&#39;s dimensions (width, height)
Definition: Kernel.h:216
int getMinX() const noexcept
Definition: Box.h:144
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
lsst::geom::Box2I shrinkBBox(lsst::geom::Box2I const &bbox) const
Given a bounding box for an image one wishes to convolve with this kernel, return the bounding box fo...
Definition: Kernel.cc:192
Reports invalid arguments.
Definition: Runtime.h:66
KernelConstPtr getKernel() const
Get the kernel (as a shared pointer to const)
Definition: Convolve.h:254
Kernels are used for convolution with MaskedImages and (eventually) Images.
Definition: Kernel.h:112
xy_locator xy_at(int x, int y) const
Return an xy_locator at the point (x, y) in the image.
Definition: ImageBase.h:443
An integer coordinate rectangle.
Definition: Box.h:54
A collection of Kernel images for special locations on a rectangular region of an image...
Definition: Convolve.h:172
virtual std::shared_ptr< Kernel > clone() const =0
Return a pointer to a deep copy of this kernel.
int getMinY() const noexcept
Definition: Box.h:145
std::ostream * os
Definition: Schema.cc:746
#define INSTANTIATE(FROMSYS, TOSYS)
Definition: Detector.cc:349