LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
Classes | Functions
lsst::afw::math::detail Namespace Reference

Classes

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...
 
struct  ConvolveWithInterpolationWorkingImages
 kernel images used by convolveRegionWithInterpolation More...
 
class  Spline
 
class  TautSpline
 
class  SmoothedSpline
 
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. More...
 
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. More...
 
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. More...
 
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. More...
 
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. More...
 
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. More...
 
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. More...
 
void declareConvolve (lsst::utils::python::WrapperCollection &wrappers)
 
void wrapConvolve (lsst::utils::python::WrapperCollection &wrappers)
 
void wrapSpline (lsst::utils::python::WrapperCollection &)
 
 PYBIND11_MODULE (_detail, mod)
 

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]convolvedImageconvolved image
[in]inImageimage to convolve
[in]kernelconvolution kernel
[in]convolutionControlconvolution control parameters

Definition at line 181 of file BasicConvolve.cc.

183  {
184  assert(!kernel.isSpatiallyVarying());
185  assertDimensionsOK(convolvedImage, inImage, kernel);
186 
187  int const mImageWidth = inImage.getWidth(); // size of input region
188  int const mImageHeight = inImage.getHeight();
189  int const cnvWidth = mImageWidth + 1 - kernel.getWidth();
190  int const cnvHeight = mImageHeight + 1 - kernel.getHeight();
191  int const cnvStartX = kernel.getCtr().getX();
192  int const cnvStartY = kernel.getCtr().getY();
193  int const inStartX = kernel.getPixel().getX();
194  int const inStartY = kernel.getPixel().getY();
195 
196  LOGL_DEBUG("TRACE2.lsst.afw.math.convolve.basicConvolve", "DeltaFunctionKernel basicConvolve");
197 
198  for (int i = 0; i < cnvHeight; ++i) {
199  typename InImageT::x_iterator inPtr = inImage.x_at(inStartX, i + inStartY);
200  for (typename OutImageT::x_iterator cnvPtr = convolvedImage.x_at(cnvStartX, i + cnvStartY),
201  cnvEnd = cnvPtr + cnvWidth;
202  cnvPtr != cnvEnd; ++cnvPtr, ++inPtr) {
203  *cnvPtr = *inPtr;
204  }
205  }
206 }
#define LOGL_DEBUG(logger, message...)
Log a debug-level message using a varargs/printf style interface.
Definition: Log.h:515

◆ 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]convolvedImageconvolved image
[in]inImageimage to convolve
[in]kernelconvolution kernel
[in]convolutionControlconvolution control parameters
Exceptions
lsst::pex::exceptions::InvalidParameterErrorif convolvedImage dimensions != inImage dimensions
lsst::pex::exceptions::InvalidParameterErrorif inImage smaller than kernel in width or height
lsst::pex::exceptions::InvalidParameterErrorif kernel width or height < 1
std::bad_allocwhen allocation of CPU memory fails

Definition at line 141 of file BasicConvolve.cc.

142  {
143  // Because convolve isn't a method of Kernel we can't always use Kernel's vtbl to dynamically
144  // dispatch the correct version of basicConvolve. The case that fails is convolving with a kernel
145  // obtained from a pointer or reference to a Kernel (base class), e.g. as used in LinearCombinationKernel.
146  if (IS_INSTANCE(kernel, math::DeltaFunctionKernel)) {
147  LOGL_DEBUG("TRACE3.lsst.afw.math.convolve.basicConvolve",
148  "generic basicConvolve: dispatch to DeltaFunctionKernel basicConvolve");
149  basicConvolve(convolvedImage, inImage, *dynamic_cast<math::DeltaFunctionKernel const*>(&kernel),
150  convolutionControl);
151  return;
152  } else if (IS_INSTANCE(kernel, math::SeparableKernel)) {
153  LOGL_DEBUG("TRACE3.lsst.afw.math.convolve.basicConvolve",
154  "generic basicConvolve: dispatch to SeparableKernel basicConvolve");
155  basicConvolve(convolvedImage, inImage, *dynamic_cast<math::SeparableKernel const*>(&kernel),
156  convolutionControl);
157  return;
158  } else if (IS_INSTANCE(kernel, math::LinearCombinationKernel) && kernel.isSpatiallyVarying()) {
159  LOGL_DEBUG(
160  "TRACE3.afw.math.convolve.basicConvolve",
161  "generic basicConvolve: dispatch to spatially varying LinearCombinationKernel basicConvolve");
162  basicConvolve(convolvedImage, inImage, *dynamic_cast<math::LinearCombinationKernel const*>(&kernel),
163  convolutionControl);
164  return;
165  }
166  // OK, use general (and slower) form
167  if (kernel.isSpatiallyVarying() && (convolutionControl.getMaxInterpolationDistance() > 1)) {
168  // use linear interpolation
169  LOGL_DEBUG("TRACE2.lsst.afw.math.convolve.basicConvolve",
170  "generic basicConvolve: using linear interpolation");
171  convolveWithInterpolation(convolvedImage, inImage, kernel, convolutionControl);
172 
173  } else {
174  // use brute force
175  LOGL_DEBUG("TRACE2.lsst.afw.math.convolve.basicConvolve", "generic basicConvolve: using brute force");
176  convolveWithBruteForce(convolvedImage, inImage, kernel, convolutionControl);
177  }
178 }
#define IS_INSTANCE(A, B)
Definition: Convolve.h:40
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.
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.
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.

◆ 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]convolvedImageconvolved image
[in]inImageimage to convolve
[in]kernelconvolution kernel
[in]convolutionControlconvolution control parameters
Exceptions
lsst::pex::exceptions::InvalidParameterErrorif convolvedImage dimensions != inImage dimensions
lsst::pex::exceptions::InvalidParameterErrorif inImage smaller than kernel in width or height
lsst::pex::exceptions::InvalidParameterErrorif kernel width or height < 1
std::bad_allocwhen allocation of CPU memory fails

Definition at line 209 of file BasicConvolve.cc.

211  {
212  if (!kernel.isSpatiallyVarying()) {
213  // use the standard algorithm for the spatially invariant case
214  LOGL_DEBUG("TRACE2.lsst.afw.math.convolve.basicConvolve",
215  "basicConvolve for LinearCombinationKernel: spatially invariant; using brute force");
216  return convolveWithBruteForce(convolvedImage, inImage, kernel, convolutionControl.getDoNormalize());
217  } else {
218  // refactor the kernel if this is reasonable and possible;
219  // then use the standard algorithm for the spatially varying case
220  std::shared_ptr<Kernel> refKernelPtr; // possibly refactored version of kernel
221  if (static_cast<int>(kernel.getNKernelParameters()) > kernel.getNSpatialParameters()) {
222  // refactoring will speed convolution, so try it
223  refKernelPtr = kernel.refactor();
224  if (!refKernelPtr) {
225  refKernelPtr = kernel.clone();
226  }
227  } else {
228  // too few basis kernels for refactoring to be worthwhile
229  refKernelPtr = kernel.clone();
230  }
231  if (convolutionControl.getMaxInterpolationDistance() > 1) {
232  LOGL_DEBUG("TRACE2.lsst.afw.math.convolve.basicConvolve",
233  "basicConvolve for LinearCombinationKernel: using interpolation");
234  return convolveWithInterpolation(convolvedImage, inImage, *refKernelPtr, convolutionControl);
235  } else {
236  LOGL_DEBUG("TRACE2.lsst.afw.math.convolve.basicConvolve",
237  "basicConvolve for LinearCombinationKernel: maxInterpolationError < 0; using brute "
238  "force");
239  return convolveWithBruteForce(convolvedImage, inImage, *refKernelPtr,
240  convolutionControl.getDoNormalize());
241  }
242  }
243 }

◆ 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]convolvedImageconvolved image
[in]inImageimage to convolve
[in]kernelconvolution kernel
[in]convolutionControlconvolution control parameters

Definition at line 246 of file BasicConvolve.cc.

247  {
248  using KernelPixel = typename math::Kernel::Pixel;
249  using KernelVector = typename std::vector<KernelPixel>;
250  using KernelIterator = KernelVector::const_iterator;
251  using InXIterator = typename InImageT::const_x_iterator;
252  using InXYLocator = typename InImageT::const_xy_locator;
253  using OutXIterator = typename OutImageT::x_iterator;
254  using OutYIterator = typename OutImageT::y_iterator;
255  using OutPixel = typename OutImageT::SinglePixel;
256 
257  assertDimensionsOK(convolvedImage, inImage, kernel);
258 
259  lsst::geom::Box2I const fullBBox = inImage.getBBox(image::LOCAL);
260  lsst::geom::Box2I const goodBBox = kernel.shrinkBBox(fullBBox);
261 
262  KernelVector kernelXVec(kernel.getWidth());
263  KernelVector kernelYVec(kernel.getHeight());
264 
265  if (kernel.isSpatiallyVarying()) {
266  LOGL_DEBUG("TRACE2.lsst.afw.math.convolve.basicConvolve",
267  "SeparableKernel basicConvolve: kernel is spatially varying");
268 
269  for (int cnvY = goodBBox.getMinY(); cnvY <= goodBBox.getMaxY(); ++cnvY) {
270  double const rowPos = inImage.indexToPosition(cnvY, image::Y);
271 
272  InXYLocator inImLoc = inImage.xy_at(0, cnvY - goodBBox.getMinY());
273  OutXIterator cnvXIter = convolvedImage.row_begin(cnvY) + goodBBox.getMinX();
274  for (int cnvX = goodBBox.getMinX(); cnvX <= goodBBox.getMaxX();
275  ++cnvX, ++inImLoc.x(), ++cnvXIter) {
276  double const colPos = inImage.indexToPosition(cnvX, image::X);
277 
278  KernelPixel kSum = kernel.computeVectors(kernelXVec, kernelYVec,
279  convolutionControl.getDoNormalize(), colPos, rowPos);
280 
281  // why does this trigger warnings? It did not in the past.
282  *cnvXIter = math::convolveAtAPoint<OutImageT, InImageT>(inImLoc, kernelXVec, kernelYVec);
283  if (convolutionControl.getDoNormalize()) {
284  *cnvXIter = *cnvXIter / kSum;
285  }
286  }
287  }
288  } else {
289  // kernel is spatially invariant
290  // The basic sequence:
291  // - For each output row:
292  // - Compute x-convolved data: a kernel height's strip of input image convolved with kernel x vector
293  // - Compute one row of output by dotting each column of x-convolved data with the kernel y vector
294  // The x-convolved data is stored in a kernel-height by good-width buffer.
295  // This is circular buffer along y (to avoid shifting pixels before setting each new row);
296  // so for each new row the kernel y vector is rotated to match the order of the x-convolved data.
297 
298  LOGL_DEBUG("TRACE2.lsst.afw.math.convolve.basicConvolve",
299  "SeparableKernel basicConvolve: kernel is spatially invariant");
300 
301  kernel.computeVectors(kernelXVec, kernelYVec, convolutionControl.getDoNormalize());
302  KernelIterator const kernelXVecBegin = kernelXVec.begin();
303  KernelIterator const kernelYVecBegin = kernelYVec.begin();
304 
305  // buffer for x-convolved data
306  OutImageT buffer(lsst::geom::Extent2I(goodBBox.getWidth(), kernel.getHeight()));
307 
308  // pre-fill x-convolved data buffer with all but one row of data
309  int yInd = 0; // during initial fill bufY = inImageY
310  int const yPrefillEnd = buffer.getHeight() - 1;
311  for (; yInd < yPrefillEnd; ++yInd) {
312  OutXIterator bufXIter = buffer.x_at(0, yInd);
313  OutXIterator const bufXEnd = buffer.x_at(goodBBox.getWidth(), yInd);
314  InXIterator inXIter = inImage.x_at(0, yInd);
315  for (; bufXIter != bufXEnd; ++bufXIter, ++inXIter) {
316  *bufXIter = kernelDotProduct<OutPixel, InXIterator, KernelIterator, KernelPixel>(
317  inXIter, kernelXVecBegin, kernel.getWidth());
318  }
319  }
320 
321  // compute output pixels using the sequence described above
322  int inY = yPrefillEnd;
323  int bufY = yPrefillEnd;
324  int cnvY = goodBBox.getMinY();
325  while (true) {
326  // fill next buffer row and compute output row
327  InXIterator inXIter = inImage.x_at(0, inY);
328  OutXIterator bufXIter = buffer.x_at(0, bufY);
329  OutXIterator cnvXIter = convolvedImage.x_at(goodBBox.getMinX(), cnvY);
330  for (int bufX = 0; bufX < goodBBox.getWidth(); ++bufX, ++cnvXIter, ++bufXIter, ++inXIter) {
331  // note: bufXIter points to the row of the buffer that is being updated,
332  // whereas bufYIter points to row 0 of the buffer
333  *bufXIter = kernelDotProduct<OutPixel, InXIterator, KernelIterator, KernelPixel>(
334  inXIter, kernelXVecBegin, kernel.getWidth());
335 
336  OutYIterator bufYIter = buffer.y_at(bufX, 0);
337  *cnvXIter = kernelDotProduct<OutPixel, OutYIterator, KernelIterator, KernelPixel>(
338  bufYIter, kernelYVecBegin, kernel.getHeight());
339  }
340 
341  // test for done now, instead of the start of the loop,
342  // to avoid an unnecessary extra rotation of the kernel Y vector
343  if (cnvY >= goodBBox.getMaxY()) break;
344 
345  // update y indices, including bufY, and rotate the kernel y vector to match
346  ++inY;
347  bufY = (bufY + 1) % kernel.getHeight();
348  ++cnvY;
349  std::rotate(kernelYVec.begin(), kernelYVec.end() - 1, kernelYVec.end());
350  }
351  }
352 }
An integer coordinate rectangle.
Definition: Box.h:55
int getMinY() const noexcept
Definition: Box.h:158
int getMinX() const noexcept
Definition: Box.h:157
int getWidth() const noexcept
Definition: Box.h:187
int getMaxX() const noexcept
Definition: Box.h:161
int getMaxY() const noexcept
Definition: Box.h:162
float Pixel
Typedefs to be used for pixel values.
Definition: common.h:37
T rotate(T... args)

◆ 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]outImageconvolved image = inImage convolved with kernel
[in]inImageinput image
[in]regionkernel image region over which to convolve
[in]workingImagesworking kernel images
Warning
: this is a low-level routine that performs no bounds checking.

Definition at line 88 of file ConvolveWithInterpolation.cc.

90  {
91  using OutLocator = typename OutImageT::xy_locator;
92  using InConstLocator = typename InImageT::const_xy_locator;
93  using KernelImage = KernelImagesForRegion::Image;
94  using KernelConstLocator = KernelImage::const_xy_locator;
95 
96  std::shared_ptr<Kernel const> kernelPtr = region.getKernel();
97  lsst::geom::Extent2I const kernelDimensions(kernelPtr->getDimensions());
98  workingImages.leftImage.assign(*region.getImage(KernelImagesForRegion::BOTTOM_LEFT));
99  workingImages.rightImage.assign(*region.getImage(KernelImagesForRegion::BOTTOM_RIGHT));
100  workingImages.kernelImage.assign(workingImages.leftImage);
101 
102  lsst::geom::Box2I const goodBBox = region.getBBox();
103  lsst::geom::Box2I const fullBBox = kernelPtr->growBBox(goodBBox);
104 
105  // top and right images are computed one beyond bbox boundary,
106  // so the distance between edge images is bbox width/height pixels
107  double xfrac = 1.0 / static_cast<double>(goodBBox.getWidth());
108  double yfrac = 1.0 / static_cast<double>(goodBBox.getHeight());
109  math::scaledPlus(workingImages.leftDeltaImage, yfrac, *region.getImage(KernelImagesForRegion::TOP_LEFT),
110  -yfrac, workingImages.leftImage);
111  math::scaledPlus(workingImages.rightDeltaImage, yfrac, *region.getImage(KernelImagesForRegion::TOP_RIGHT),
112  -yfrac, workingImages.rightImage);
113 
114  KernelConstLocator const kernelLocator = workingImages.kernelImage.xy_at(0, 0);
115 
116  // The loop is a bit odd for efficiency: the initial value of workingImages.kernelImage
117  // and related kernel images are set when they are allocated,
118  // so they are not computed in the loop until after the convolution; to save cpu cycles
119  // they are not computed at all for the last iteration.
120  InConstLocator inLocator = inImage.xy_at(fullBBox.getMinX(), fullBBox.getMinY());
121  OutLocator outLocator = outImage.xy_at(goodBBox.getMinX(), goodBBox.getMinY());
122  for (int j = 0;;) {
123  auto inLocatorInitialPosition = inLocator;
124  auto outLocatorInitialPosition = outLocator;
125  math::scaledPlus(workingImages.deltaImage, xfrac, workingImages.rightImage, -xfrac,
126  workingImages.leftImage);
127  for (int i = 0;;) {
128  *outLocator = math::convolveAtAPoint<OutImageT, InImageT>(
129  inLocator, kernelLocator, kernelDimensions.getX(), kernelDimensions.getY());
130  ++outLocator.x();
131  ++inLocator.x();
132  ++i;
133  if (i >= goodBBox.getWidth()) {
134  break;
135  }
136  workingImages.kernelImage += workingImages.deltaImage;
137  }
138 
139  ++j;
140  if (j >= goodBBox.getHeight()) {
141  break;
142  }
143  workingImages.leftImage += workingImages.leftDeltaImage;
144  workingImages.rightImage += workingImages.rightDeltaImage;
145  workingImages.kernelImage.assign(workingImages.leftImage);
146 
147  // Workaround for DM-5822
148  // Boost GIL locator won't decrement in x-dimension for some strange and still
149  // not understood reason. So instead store position at start of previous row,
150  // reset and move down.
151  inLocator = inLocatorInitialPosition;
152  outLocator = outLocatorInitialPosition;
153  ++inLocator.y();
154  ++outLocator.y();
155  }
156 }
int getHeight() const noexcept
Definition: Box.h:188
void scaledPlus(OutImageT &outImage, double c1, InImageT const &inImage1, double c2, InImageT const &inImage2)
Compute the scaled sum of two images.

◆ 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]convolvedImageconvolved image
[in]inImageimage to convolve
[in]kernelconvolution kernel
[in]convolutionControlconvolution control parameters
Exceptions
lsst::pex::exceptions::InvalidParameterErrorif convolvedImage dimensions != inImage dimensions
lsst::pex::exceptions::InvalidParameterErrorif inImage smaller than kernel in width or height
lsst::pex::exceptions::InvalidParameterErrorif kernel width or height < 1
std::bad_allocwhen allocation of CPU memory fails
Warning
Low-level convolution function that does not set edge pixels.

Definition at line 355 of file BasicConvolve.cc.

356  {
357  bool doNormalize = convolutionControl.getDoNormalize();
358 
359  using KernelPixel = typename math::Kernel::Pixel;
360  using KernelImage = image::Image<KernelPixel>;
361 
362  using KernelXIterator = typename KernelImage::const_x_iterator;
363  using KernelXYLocator = typename KernelImage::const_xy_locator;
364  using InXIterator = typename InImageT::const_x_iterator;
365  using InXYLocator = typename InImageT::const_xy_locator;
366  using OutXIterator = typename OutImageT::x_iterator;
367  using OutPixel = typename OutImageT::SinglePixel;
368 
369  assertDimensionsOK(convolvedImage, inImage, kernel);
370 
371  int const inImageWidth = inImage.getWidth();
372  int const inImageHeight = inImage.getHeight();
373  int const kWidth = kernel.getWidth();
374  int const kHeight = kernel.getHeight();
375  int const cnvWidth = inImageWidth + 1 - kernel.getWidth();
376  int const cnvHeight = inImageHeight + 1 - kernel.getHeight();
377  int const cnvStartX = kernel.getCtr().getX();
378  int const cnvStartY = kernel.getCtr().getY();
379  int const cnvEndX = cnvStartX + cnvWidth; // end index + 1
380  int const cnvEndY = cnvStartY + cnvHeight; // end index + 1
381 
382  KernelImage kernelImage(kernel.getDimensions());
383  KernelXYLocator const kernelLoc = kernelImage.xy_at(0, 0);
384 
385  if (kernel.isSpatiallyVarying()) {
386  LOGL_DEBUG("TRACE4.lsst.afw.math.convolve.convolveWithBruteForce",
387  "convolveWithBruteForce: kernel is spatially varying");
388 
389  for (int cnvY = cnvStartY; cnvY != cnvEndY; ++cnvY) {
390  double const rowPos = inImage.indexToPosition(cnvY, image::Y);
391 
392  InXYLocator inImLoc = inImage.xy_at(0, cnvY - cnvStartY);
393  OutXIterator cnvXIter = convolvedImage.x_at(cnvStartX, cnvY);
394  for (int cnvX = cnvStartX; cnvX != cnvEndX; ++cnvX, ++inImLoc.x(), ++cnvXIter) {
395  double const colPos = inImage.indexToPosition(cnvX, image::X);
396 
397  KernelPixel kSum = kernel.computeImage(kernelImage, false, colPos, rowPos);
398  *cnvXIter = math::convolveAtAPoint<OutImageT, InImageT>(inImLoc, kernelLoc, kWidth, kHeight);
399  if (doNormalize) {
400  *cnvXIter = *cnvXIter / kSum;
401  }
402  }
403  }
404  } else {
405  LOGL_DEBUG("TRACE4.lsst.afw.math.convolve.convolveWithBruteForce",
406  "convolveWithBruteForce: kernel is spatially invariant");
407 
408  (void)kernel.computeImage(kernelImage, doNormalize);
409 
410  for (int inStartY = 0, cnvY = cnvStartY; inStartY < cnvHeight; ++inStartY, ++cnvY) {
411  KernelXIterator kernelXIter = kernelImage.x_at(0, 0);
412  InXIterator inXIter = inImage.x_at(0, inStartY);
413  OutXIterator cnvXIter = convolvedImage.x_at(cnvStartX, cnvY);
414  for (int x = 0; x < cnvWidth; ++x, ++cnvXIter, ++inXIter) {
415  *cnvXIter = kernelDotProduct<OutPixel, InXIterator, KernelXIterator, KernelPixel>(
416  inXIter, kernelXIter, kWidth);
417  }
418  for (int kernelY = 1, inY = inStartY + 1; kernelY < kHeight; ++inY, ++kernelY) {
419  KernelXIterator kernelXIter = kernelImage.x_at(0, kernelY);
420  InXIterator inXIter = inImage.x_at(0, inY);
421  OutXIterator cnvXIter = convolvedImage.x_at(cnvStartX, cnvY);
422  for (int x = 0; x < cnvWidth; ++x, ++cnvXIter, ++inXIter) {
423  *cnvXIter += kernelDotProduct<OutPixel, InXIterator, KernelXIterator, KernelPixel>(
424  inXIter, kernelXIter, kWidth);
425  }
426  }
427  }
428  }
429 }
double x
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:51

◆ 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]outImageconvolved image = inImage convolved with kernel
[in]inImageinput image
[in]kernelconvolution kernel
[in]convolutionControlconvolution control parameters
Exceptions
lsst::pex::exceptions::InvalidParameterErrorif outImage is not the same size as inImage

Definition at line 45 of file ConvolveWithInterpolation.cc.

46  {
47  if (outImage.getDimensions() != inImage.getDimensions()) {
49  os << "outImage dimensions = ( " << outImage.getWidth() << ", " << outImage.getHeight() << ") != ("
50  << inImage.getWidth() << ", " << inImage.getHeight() << ") = inImage dimensions";
52  }
53 
54  // compute region covering good area of output image
56  lsst::geom::Point2I(0, 0), lsst::geom::Extent2I(outImage.getWidth(), outImage.getHeight()));
57  lsst::geom::Box2I goodBBox = kernel.shrinkBBox(fullBBox);
58  KernelImagesForRegion goodRegion(KernelImagesForRegion(kernel.clone(), goodBBox, inImage.getXY0(),
59  convolutionControl.getDoNormalize()));
60  LOGL_DEBUG("TRACE5.lsst.afw.math.convolve.convolveWithInterpolation",
61  "convolveWithInterpolation: full bbox minimum=(%d, %d), extent=(%d, %d)", fullBBox.getMinX(),
62  fullBBox.getMinY(), fullBBox.getWidth(), fullBBox.getHeight());
63  LOGL_DEBUG("TRACE5.lsst.afw.math.convolve.convolveWithInterpolation",
64  "convolveWithInterpolation: goodRegion bbox minimum=(%d, %d), extent=(%d, %d)",
65  goodRegion.getBBox().getMinX(), goodRegion.getBBox().getMinY(),
66  goodRegion.getBBox().getWidth(), goodRegion.getBBox().getHeight());
67 
68  // divide good region into subregions small enough to interpolate over
69  int nx = 1 + (goodBBox.getWidth() / convolutionControl.getMaxInterpolationDistance());
70  int ny = 1 + (goodBBox.getHeight() / convolutionControl.getMaxInterpolationDistance());
71  LOGL_DEBUG("TRACE3.lsst.afw.math.convolve.convolveWithInterpolation",
72  "convolveWithInterpolation: divide into %d x %d subregions", nx, ny);
73 
74  ConvolveWithInterpolationWorkingImages workingImages(kernel.getDimensions());
75  RowOfKernelImagesForRegion regionRow(nx, ny);
76  while (goodRegion.computeNextRow(regionRow)) {
77  for (auto const &rgnIter : regionRow) {
78  LOGL_DEBUG("TRACE5.lsst.afw.math.convolve.convolveWithInterpolation",
79  "convolveWithInterpolation: bbox minimum=(%d, %d), extent=(%d, %d)",
80  rgnIter->getBBox().getMinX(), rgnIter->getBBox().getMinY(),
81  rgnIter->getBBox().getWidth(), rgnIter->getBBox().getHeight());
82  convolveRegionWithInterpolation(outImage, inImage, *rgnIter, workingImages);
83  }
84  }
85 }
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
std::ostream * os
Definition: Schema.cc:557
Reports invalid arguments.
Definition: Runtime.h:66
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.

◆ declareConvolve()

void lsst::afw::math::detail::declareConvolve ( lsst::utils::python::WrapperCollection &  wrappers)

Definition at line 69 of file _convolve.cc.

69  {
70  using PyClass = py::class_<KernelImagesForRegion, std::shared_ptr<KernelImagesForRegion>>;
71  auto clsKernelImagesForRegion =
72  wrappers.wrapType(PyClass(wrappers.module, "KernelImagesForRegion"), [](auto &mod, auto &cls) {
73  cls.def(py::init<KernelImagesForRegion::KernelConstPtr, lsst::geom::Box2I const &,
74  lsst::geom::Point2I const &, bool>(),
75  "kernelPtr"_a, "bbox"_a, "xy0"_a, "doNormalize"_a);
76  cls.def(py::init<KernelImagesForRegion::KernelConstPtr, lsst::geom::Box2I const &,
77  lsst::geom::Point2I const &, bool, KernelImagesForRegion::ImagePtr,
78  KernelImagesForRegion::ImagePtr, KernelImagesForRegion::ImagePtr,
79  KernelImagesForRegion::ImagePtr>(),
80  "kernelPtr"_a, "bbox"_a, "xy0"_a, "doNormalize"_a, "bottomLeftImagePtr"_a,
81  "bottomRightImagePtr"_a, "topLeftImagePtr"_a, "topRightImagePtr"_a);
82 
83  cls.def("getBBox", &KernelImagesForRegion::getBBox);
84  cls.def("getXY0", &KernelImagesForRegion::getXY0);
85  cls.def("getDoNormalize", &KernelImagesForRegion::getDoNormalize);
86  cls.def("getImage", &KernelImagesForRegion::getImage);
87  cls.def("getKernel", &KernelImagesForRegion::getKernel);
88  cls.def("getPixelIndex", &KernelImagesForRegion::getPixelIndex);
89  cls.def("computeNextRow", &KernelImagesForRegion::computeNextRow);
90  cls.def_static("getMinInterpolationSize", KernelImagesForRegion::getMinInterpolationSize);
91  });
92 
93  wrappers.wrapType(py::enum_<KernelImagesForRegion::Location>(clsKernelImagesForRegion, "Location"),
94  [](auto &mod, auto &enm) {
95  enm.value("BOTTOM_LEFT", KernelImagesForRegion::Location::BOTTOM_LEFT);
96  enm.value("BOTTOM_RIGHT", KernelImagesForRegion::Location::BOTTOM_RIGHT);
97  enm.value("TOP_LEFT", KernelImagesForRegion::Location::TOP_LEFT);
98  enm.value("TOP_RIGHT", KernelImagesForRegion::Location::TOP_RIGHT);
99  enm.export_values();
100  });
101 
102  wrappers.wrapType(py::class_<RowOfKernelImagesForRegion, std::shared_ptr<RowOfKernelImagesForRegion>>(
103  wrappers.module, "RowOfKernelImagesForRegion"),
104  [](auto &mod, auto &cls) {
105  cls.def(py::init<int, int>(), "nx"_a, "ny"_a);
106 
107  cls.def("front", &RowOfKernelImagesForRegion::front);
108  cls.def("back", &RowOfKernelImagesForRegion::back);
109  cls.def("getNX", &RowOfKernelImagesForRegion::getNX);
110  cls.def("getNY", &RowOfKernelImagesForRegion::getNY);
111  cls.def("getYInd", &RowOfKernelImagesForRegion::getYInd);
112  cls.def("getRegion", &RowOfKernelImagesForRegion::getRegion);
113  cls.def("hasData", &RowOfKernelImagesForRegion::hasData);
114  cls.def("isLastRow", &RowOfKernelImagesForRegion::isLastRow);
115  cls.def("incrYInd", &RowOfKernelImagesForRegion::incrYInd);
116  });
117 }
py::class_< PixelAreaBoundedField, std::shared_ptr< PixelAreaBoundedField >, BoundedField > PyClass

◆ PYBIND11_MODULE()

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

Definition at line 11 of file _detail.cc.

11  {
12  lsst::utils::python::WrapperCollection wrappers(mod, "lsst.afw.math.detail");
13  wrapConvolve(wrappers);
14  wrapSpline(wrappers);
15  wrappers.finish();
16 }
void wrapSpline(lsst::utils::python::WrapperCollection &)
Definition: _spline.cc:38
void wrapConvolve(lsst::utils::python::WrapperCollection &wrappers)
Definition: _convolve.cc:118

◆ wrapConvolve()

void lsst::afw::math::detail::wrapConvolve ( lsst::utils::python::WrapperCollection &  wrappers)

Definition at line 118 of file _convolve.cc.

118  {
119  wrappers.addSignatureDependency("lsst.afw.image");
120  wrappers.addSignatureDependency("lsst.afw.math");
121  declareConvolve(wrappers);
122  declareAll<double, double>(wrappers);
123  declareAll<double, float>(wrappers);
124  declareAll<double, int>(wrappers);
125  declareAll<double, std::uint16_t>(wrappers);
126  declareAll<float, float>(wrappers);
127  declareAll<float, int>(wrappers);
128  declareAll<float, std::uint16_t>(wrappers);
129  declareAll<int, int>(wrappers);
130  declareAll<std::uint16_t, std::uint16_t>(wrappers);
131 }
void declareConvolve(lsst::utils::python::WrapperCollection &wrappers)
Definition: _convolve.cc:69

◆ wrapSpline()

void lsst::afw::math::detail::wrapSpline ( lsst::utils::python::WrapperCollection &  wrappers)

Definition at line 38 of file _spline.cc.

38  {
39  /* Module level */
40  wrappers.wrapType(py::class_<Spline>(wrappers.module, "Spline"), [](auto &mod, auto &cls) {
41  cls.def("interpolate", &Spline::interpolate);
42  cls.def("derivative", &Spline::derivative);
43  });
44  auto clsTautSpline = wrappers.wrapType(
45  py::class_<TautSpline, Spline>(wrappers.module, "TautSpline"), [](auto &mod, auto &cls) {
46  cls.def(py::init<std::vector<double> const &, std::vector<double> const &, double const,
47  TautSpline::Symmetry>(),
48  "x"_a, "y"_a, "gamma"_a = 0, "type"_a = lsst::afw::math::detail::TautSpline::Unknown);
49  cls.def("roots", &TautSpline::roots);
50  });
51  /* Member types and enums */
52  wrappers.wrapType(py::enum_<TautSpline::Symmetry>(clsTautSpline, "Symmetry"), [](auto &mod, auto &enm) {
53  enm.value("Unknown", TautSpline::Symmetry::Unknown);
54  enm.value("Odd", TautSpline::Symmetry::Odd);
55  enm.value("Even", TautSpline::Symmetry::Even);
56  enm.export_values();
57  });
58 }