LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
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 
34 #include <algorithm>
35 #include <cmath>
36 #include <sstream>
37 #include <vector>
38 #include <iostream>
39 
40 #include "boost/cstdint.hpp"
41 
42 #include "lsst/pex/exceptions.h"
43 #include "lsst/pex/logging/Trace.h"
45 #include "lsst/afw/math/Kernel.h"
46 #include "lsst/afw/geom.h"
48 
49 namespace pexExcept = lsst::pex::exceptions;
50 namespace pexLog = lsst::pex::logging;
51 namespace afwGeom = lsst::afw::geom;
52 namespace afwImage = lsst::afw::image;
53 namespace afwMath = lsst::afw::math;
54 namespace mathDetail = lsst::afw::math::detail;
55 
70 template <typename OutImageT, typename InImageT>
72  OutImageT &outImage,
73  InImageT const &inImage,
74  lsst::afw::math::Kernel const &kernel,
75  lsst::afw::math::ConvolutionControl const &convolutionControl)
76 {
77  if (outImage.getDimensions() != inImage.getDimensions()) {
78  std::ostringstream os;
79  os << "outImage dimensions = ( "
80  << outImage.getWidth() << ", " << outImage.getHeight()
81  << ") != (" << inImage.getWidth() << ", " << inImage.getHeight() << ") = inImage dimensions";
82  throw LSST_EXCEPT(pexExcept::InvalidParameterError, os.str());
83  }
84 
85  // compute region covering good area of output image
86  afwGeom::Box2I fullBBox = afwGeom::Box2I(
87  afwGeom::Point2I(0, 0),
88  afwGeom::Extent2I(outImage.getWidth(), outImage.getHeight()));
89  afwGeom::Box2I goodBBox = kernel.shrinkBBox(fullBBox);
91  kernel.clone(),
92  goodBBox,
93  inImage.getXY0(),
94  convolutionControl.getDoNormalize()));
95  pexLog::TTrace<6>("lsst.afw.math.convolve",
96  "convolveWithInterpolation: full bbox minimum=(%d, %d), extent=(%d, %d)",
97  fullBBox.getMinX(), fullBBox.getMinY(),
98  fullBBox.getWidth(), fullBBox.getHeight());
99  pexLog::TTrace<6>("lsst.afw.math.convolve",
100  "convolveWithInterpolation: goodRegion bbox minimum=(%d, %d), extent=(%d, %d)",
101  goodRegion.getBBox().getMinX(), goodRegion.getBBox().getMinY(),
102  goodRegion.getBBox().getWidth(), goodRegion.getBBox().getHeight());
103 
104  // divide good region into subregions small enough to interpolate over
105  int nx = 1 + (goodBBox.getWidth() / convolutionControl.getMaxInterpolationDistance());
106  int ny = 1 + (goodBBox.getHeight() / convolutionControl.getMaxInterpolationDistance());
107  pexLog::TTrace<4>("lsst.afw.math.convolve",
108  "convolveWithInterpolation: divide into %d x %d subregions", nx, ny);
109 
111  RowOfKernelImagesForRegion regionRow(nx, ny);
112  while (goodRegion.computeNextRow(regionRow)) {
113  for (RowOfKernelImagesForRegion::ConstIterator rgnIter = regionRow.begin(), rgnEnd = regionRow.end();
114  rgnIter != rgnEnd; ++rgnIter) {
115  pexLog::TTrace<6>("lsst.afw.math.convolve",
116  "convolveWithInterpolation: bbox minimum=(%d, %d), extent=(%d, %d)",
117  (*rgnIter)->getBBox().getMinX(), (*rgnIter)->getBBox().getMinY(),
118  (*rgnIter)->getBBox().getWidth(), (*rgnIter)->getBBox().getHeight());
119  convolveRegionWithInterpolation(outImage, inImage, **rgnIter, workingImages);
120  }
121  }
122 }
123 
131 template <typename OutImageT, typename InImageT>
133  OutImageT &outImage,
134  InImageT const &inImage,
135  KernelImagesForRegion const &region,
137 {
138  typedef typename OutImageT::xy_locator OutLocator;
139  typedef typename InImageT::const_xy_locator InConstLocator;
140  typedef KernelImagesForRegion::Image KernelImage;
141  typedef KernelImage::const_xy_locator KernelConstLocator;
142 
143  CONST_PTR(afwMath::Kernel) kernelPtr = region.getKernel();
144  geom::Extent2I const kernelDimensions(kernelPtr->getDimensions());
145  workingImages.leftImage <<= *region.getImage(KernelImagesForRegion::BOTTOM_LEFT);
146  workingImages.rightImage <<= *region.getImage(KernelImagesForRegion::BOTTOM_RIGHT);
147  workingImages.kernelImage <<= workingImages.leftImage;
148 
149  afwGeom::Box2I const goodBBox = region.getBBox();
150  afwGeom::Box2I const fullBBox = kernelPtr->growBBox(goodBBox);
151 
152  // top and right images are computed one beyond bbox boundary,
153  // so the distance between edge images is bbox width/height pixels
154  double xfrac = 1.0 / static_cast<double>(goodBBox.getWidth());
155  double yfrac = 1.0 / static_cast<double>(goodBBox.getHeight());
156  afwMath::scaledPlus(workingImages.leftDeltaImage,
158  -yfrac, workingImages.leftImage);
159  afwMath::scaledPlus(workingImages.rightDeltaImage,
161  -yfrac, workingImages.rightImage);
162 
163  KernelConstLocator const kernelLocator = workingImages.kernelImage.xy_at(0, 0);
164 
165  // The loop is a bit odd for efficiency: the initial value of workingImages.kernelImage
166  // and related kernel images are set when they are allocated,
167  // so they are not computed in the loop until after the convolution; to save cpu cycles
168  // they are not computed at all for the last iteration.
169  InConstLocator inLocator = inImage.xy_at(fullBBox.getMinX(), fullBBox.getMinY());
170  OutLocator outLocator = outImage.xy_at(goodBBox.getMinX(), goodBBox.getMinY());
171  for (int j = 0; ; ) {
173  workingImages.deltaImage, xfrac, workingImages.rightImage, -xfrac, workingImages.leftImage);
174  for (int i = 0; ; ) {
175  *outLocator = afwMath::convolveAtAPoint<OutImageT, InImageT>(
176  inLocator, kernelLocator, kernelDimensions.getX(), kernelDimensions.getY());
177  ++outLocator.x();
178  ++inLocator.x();
179  ++i;
180  if (i >= goodBBox.getWidth()) {
181  break;
182  }
183  workingImages.kernelImage += workingImages.deltaImage;
184  }
185 
186  ++j;
187  if (j >= goodBBox.getHeight()) {
188  break;
189  }
190  workingImages.leftImage += workingImages.leftDeltaImage;
191  workingImages.rightImage += workingImages.rightDeltaImage;
192  workingImages.kernelImage <<= workingImages.leftImage;
193  inLocator += lsst::afw::image::detail::difference_type(-goodBBox.getWidth(), 1);
194  outLocator += lsst::afw::image::detail::difference_type(-goodBBox.getWidth(), 1);
195  }
196 }
197 
198 /*
199  * Explicit instantiation
200  */
202 #define IMAGE(PIXTYPE) afwImage::Image<PIXTYPE>
203 #define MASKEDIMAGE(PIXTYPE) afwImage::MaskedImage<PIXTYPE, afwImage::MaskPixel, afwImage::VariancePixel>
204 #define NL /* */
205 // Instantiate Image or MaskedImage versions
206 #define INSTANTIATE_IM_OR_MI(IMGMACRO, OUTPIXTYPE, INPIXTYPE) \
207  template void mathDetail::convolveWithInterpolation( \
208  IMGMACRO(OUTPIXTYPE)&, IMGMACRO(INPIXTYPE) const&, afwMath::Kernel const&, \
209  afwMath::ConvolutionControl const&); NL \
210  template void mathDetail::convolveRegionWithInterpolation( \
211  IMGMACRO(OUTPIXTYPE)&, IMGMACRO(INPIXTYPE) const&, KernelImagesForRegion const&, \
212  ConvolveWithInterpolationWorkingImages&);
213 // Instantiate both Image and MaskedImage versions
214 #define INSTANTIATE(OUTPIXTYPE, INPIXTYPE) \
215  INSTANTIATE_IM_OR_MI(IMAGE, OUTPIXTYPE, INPIXTYPE) \
216  INSTANTIATE_IM_OR_MI(MASKEDIMAGE, OUTPIXTYPE, INPIXTYPE)
217 
218 INSTANTIATE(double, double)
219 INSTANTIATE(double, float)
220 INSTANTIATE(double, int)
221 INSTANTIATE(double, boost::uint16_t)
222 INSTANTIATE(float, float)
223 INSTANTIATE(float, int)
224 INSTANTIATE(float, boost::uint16_t)
225 INSTANTIATE(int, int)
226 INSTANTIATE(boost::uint16_t, boost::uint16_t)
Convolution support.
void scaledPlus(OutImageT &outImage, double c1, InImageT const &inImage1, double c2, InImageT const &inImage2)
An include file to include the header files for lsst::afw::geom.
Declare the Kernel class and subclasses.
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:58
#define CONST_PTR(...)
Definition: base.h:47
definition of the Trace messaging facilities
lsst::afw::geom::Box2I getBBox() const
Definition: Convolve.h:158
int getMinY() const
Definition: Box.h:125
ImagePtr getImage(Location location) const
An integer coordinate rectangle.
Definition: Box.h:53
kernel images used by convolveRegionWithInterpolation
Definition: Convolve.h:266
table::Key< table::Array< Kernel::Pixel > > image
Definition: FixedKernel.cc:117
virtual boost::shared_ptr< Kernel > clone() const =0
Return a pointer to a deep copy of this kernel.
int getHeight() const
Definition: Box.h:155
geom::Extent2I const getDimensions() const
Return the Kernel&#39;s dimensions (width, height)
Definition: Kernel.h:226
A row of KernelImagesForRegion.
Definition: Convolve.h:205
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.
int getMinX() const
Definition: Box.h:124
bool computeNextRow(RowOfKernelImagesForRegion &regionRow) const
Compute next row of subregions.
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
#define INSTANTIATE(T)
xy_locator xy_at(int x, int y) const
Definition: Image.h:351
Implementation of the Class MaskedImage.
RegionList::const_iterator ConstIterator
Definition: Convolve.h:209
lsst::afw::geom::Box2I shrinkBBox(lsst::afw::geom::Box2I const &bbox) const
Definition: Kernel.cc:207
Kernels are used for convolution with MaskedImages and (eventually) Images.
Definition: Kernel.h:134
int getWidth() const
Definition: Box.h:154
Include files required for standard LSST Exception handling.