LSST Applications  21.0.0+04719a4bac,21.0.0-1-ga51b5d4+f5e6047307,21.0.0-11-g2b59f77+a9c1acf22d,21.0.0-11-ga42c5b2+86977b0b17,21.0.0-12-gf4ce030+76814010d2,21.0.0-13-g1721dae+760e7a6536,21.0.0-13-g3a573fe+768d78a30a,21.0.0-15-g5a7caf0+f21cbc5713,21.0.0-16-g0fb55c1+b60e2d390c,21.0.0-19-g4cded4ca+71a93a33c0,21.0.0-2-g103fe59+bb20972958,21.0.0-2-g45278ab+04719a4bac,21.0.0-2-g5242d73+3ad5d60fb1,21.0.0-2-g7f82c8f+8babb168e8,21.0.0-2-g8f08a60+06509c8b61,21.0.0-2-g8faa9b5+616205b9df,21.0.0-2-ga326454+8babb168e8,21.0.0-2-gde069b7+5e4aea9c2f,21.0.0-2-gecfae73+1d3a86e577,21.0.0-2-gfc62afb+3ad5d60fb1,21.0.0-25-g1d57be3cd+e73869a214,21.0.0-3-g357aad2+ed88757d29,21.0.0-3-g4a4ce7f+3ad5d60fb1,21.0.0-3-g4be5c26+3ad5d60fb1,21.0.0-3-g65f322c+e0b24896a3,21.0.0-3-g7d9da8d+616205b9df,21.0.0-3-ge02ed75+a9c1acf22d,21.0.0-4-g591bb35+a9c1acf22d,21.0.0-4-g65b4814+b60e2d390c,21.0.0-4-gccdca77+0de219a2bc,21.0.0-4-ge8a399c+6c55c39e83,21.0.0-5-gd00fb1e+05fce91b99,21.0.0-6-gc675373+3ad5d60fb1,21.0.0-64-g1122c245+4fb2b8f86e,21.0.0-7-g04766d7+cd19d05db2,21.0.0-7-gdf92d54+04719a4bac,21.0.0-8-g5674e7b+d1bd76f71f,master-gac4afde19b+a9c1acf22d,w.2021.13
LSST Data Management Base Package
warpExposure.cc
Go to the documentation of this file.
1 // -*- LSST-C++ -*- // fixed format comment for emacs
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  * Support for warping an %image to a new Wcs.
26  */
27 
28 #include <cassert>
29 #include <cmath>
30 #include <cstdint>
31 #include <limits>
32 #include <memory>
33 #include <sstream>
34 #include <string>
35 #include <vector>
36 #include <utility>
37 #include <ctime>
38 
39 #include <memory>
40 #include "boost/pointer_cast.hpp"
41 #include "boost/regex.hpp"
42 #include "astshim.h"
43 
44 #include "lsst/log/Log.h"
45 #include "lsst/pex/exceptions.h"
46 #include "lsst/geom.h"
48 #include "lsst/afw/geom.h"
49 #include "lsst/afw/math/Kernel.h"
52 
54 
55 using std::swap;
56 
57 namespace lsst {
58 namespace afw {
59 namespace math {
60 
61 //
62 // A helper function for the warping kernels which provides error-checking:
63 // the warping kernels are designed to work in two cases
64 // 0 < x < 1 and ctrX=(size-1)/2
65 // -1 < x < 0 and ctrX=(size+1)/2
66 // (and analogously for y). Note that to get the second case, Kernel::setCtr(1, y) must be
67 // called before calling Kernel::setKernelParameter(). [see afw::math::offsetImage() for
68 // an example]
69 //
70 // FIXME eventually the 3 warping kernels will inherit from a common base class WarpingKernel
71 // and this routine can be eliminated by putting the code in WarpingKernel::setKernelParameter()
72 //
73 static inline void checkWarpingKernelParameter(const SeparableKernel *p, unsigned int ind, double value) {
74  if (ind > 1) {
76  "bad ind argument in WarpingKernel::setKernelParameter()");
77  }
78  int ctr = p->getCtr()[ind];
79  int size = p->getDimensions()[ind];
80 
81  if (ctr == (size - 1) / 2) {
82  if (value < -1e-6 || value > 1 + 1e-6) {
84  "bad coordinate in WarpingKernel::setKernelParameter()");
85  }
86  } else if (ctr == (size + 1) / 2) {
87  if (value < -1 - 1e-6 || value > 1e-6) {
89  "bad coordinate in WarpingKernel::setKernelParameter()");
90  }
91  } else {
93  "bad ctr value in WarpingKernel::setKernelParameter()");
94  }
95 }
96 
99 }
100 
101 int LanczosWarpingKernel::getOrder() const { return this->getWidth() / 2; }
102 
103 void LanczosWarpingKernel::setKernelParameter(unsigned int ind, double value) const {
104  checkWarpingKernelParameter(this, ind, value);
106 }
107 
110 }
111 
113  //
114  // this->_params[0] = value of x where we want to interpolate the function
115  // x = integer value of x where we evaluate the function in the interpolation
116  //
117  // The following weird-looking expression has no if/else statements, is roundoff-tolerant,
118  // and works in the following two cases:
119  // 0 < this->_params[0] < 1, x \in {0,1}
120  // -1 < this->_params[0] < 0, x \in {-1,0}
121  //
122  // The checks in BilinearWarpingKernel::setKernelParameter() ensure that one of these
123  // conditions is satisfied
124  //
125  return 0.5 + (1.0 - (2.0 * fabs(this->_params[0]))) * (0.5 - fabs(x));
126 }
127 
128 void BilinearWarpingKernel::setKernelParameter(unsigned int ind, double value) const {
129  checkWarpingKernelParameter(this, ind, value);
131 }
132 
135  os << "_BilinearFunction1: ";
136  os << Function1<Kernel::Pixel>::toString(prefix);
137  return os.str();
138 }
139 
141  return std::make_shared<NearestWarpingKernel>();
142 }
143 
145  // this expression is faster than using conditionals, but offers no sanity checking
146  return static_cast<double>((fabs(this->_params[0]) < 0.5) == (fabs(x) < 0.5));
147 }
148 
149 void NearestWarpingKernel::setKernelParameter(unsigned int ind, double value) const {
150  checkWarpingKernelParameter(this, ind, value);
152 }
153 
156  os << "_NearestFunction1: ";
157  os << Function1<Kernel::Pixel>::toString(prefix);
158  return os.str();
159 }
160 
162  typedef std::shared_ptr<SeparableKernel> KernelPtr;
163  boost::cmatch matches;
164  static const boost::regex LanczosRE("lanczos(\\d+)");
165  if (name == "bilinear") {
166  return KernelPtr(new BilinearWarpingKernel());
167  } else if (boost::regex_match(name.c_str(), matches, LanczosRE)) {
168  std::string orderStr(matches[1].first, matches[1].second);
169  int order;
170  std::istringstream(orderStr) >> order;
171  return KernelPtr(new LanczosWarpingKernel(order));
172  } else if (name == "nearest") {
173  return KernelPtr(new NearestWarpingKernel());
174  } else {
175  throw LSST_EXCEPT(pexExcept::InvalidParameterError, "unknown warping kernel name: \"" + name + "\"");
176  }
177 }
178 
180  if (_warpingKernelPtr->getCacheSize() != _cacheSize) {
181  _warpingKernelPtr->computeCache(_cacheSize);
182  }
183  return _warpingKernelPtr;
184 };
185 
188  setWarpingKernel(*warpingKernelPtr);
189 }
190 
192  if (_maskWarpingKernelPtr) {
193  _testWarpingKernels(warpingKernel, *_maskWarpingKernelPtr);
194  }
195  std::shared_ptr<SeparableKernel> warpingKernelPtr(
196  std::static_pointer_cast<SeparableKernel>(warpingKernel.clone()));
197  _warpingKernelPtr = warpingKernelPtr;
198 }
199 
201  if (_maskWarpingKernelPtr) { // lazily update kernel cache
202  if (_maskWarpingKernelPtr->getCacheSize() != _cacheSize) {
203  _maskWarpingKernelPtr->computeCache(_cacheSize);
204  }
205  }
206  return _maskWarpingKernelPtr;
207 }
208 
209 void WarpingControl::setMaskWarpingKernelName(std::string const &maskWarpingKernelName) {
210  if (!maskWarpingKernelName.empty()) {
211  std::shared_ptr<SeparableKernel> maskWarpingKernelPtr(makeWarpingKernel(maskWarpingKernelName));
212  setMaskWarpingKernel(*maskWarpingKernelPtr);
213  } else {
214  _maskWarpingKernelPtr.reset();
215  }
216 }
217 
219  _testWarpingKernels(*_warpingKernelPtr, maskWarpingKernel);
220  _maskWarpingKernelPtr = std::static_pointer_cast<SeparableKernel>(maskWarpingKernel.clone());
221 }
222 
223 void WarpingControl::_testWarpingKernels(SeparableKernel const &warpingKernel,
224  SeparableKernel const &maskWarpingKernel) const {
225  lsst::geom::Box2I kernelBBox =
227  warpingKernel.getDimensions());
228  lsst::geom::Box2I maskKernelBBox =
230  maskWarpingKernel.getDimensions());
231  if (!kernelBBox.contains(maskKernelBBox)) {
233  "warping kernel is smaller than mask warping kernel");
234  }
235 }
236 
237 template <typename DestExposureT, typename SrcExposureT>
238 int warpExposure(DestExposureT &destExposure, SrcExposureT const &srcExposure, WarpingControl const &control,
239  typename DestExposureT::MaskedImageT::SinglePixel padValue) {
240  if (!destExposure.hasWcs()) {
241  throw LSST_EXCEPT(pexExcept::InvalidParameterError, "destExposure has no Wcs");
242  }
243  if (!srcExposure.hasWcs()) {
244  throw LSST_EXCEPT(pexExcept::InvalidParameterError, "srcExposure has no Wcs");
245  }
246  typename DestExposureT::MaskedImageT mi = destExposure.getMaskedImage();
247  destExposure.setPhotoCalib(srcExposure.getPhotoCalib());
248  destExposure.setFilterLabel(srcExposure.getFilterLabel());
249  destExposure.getInfo()->setVisitInfo(srcExposure.getInfo()->getVisitInfo());
250  return warpImage(mi, *destExposure.getWcs(), srcExposure.getMaskedImage(), *srcExposure.getWcs(), control,
251  padValue);
252 }
253 
254 namespace {
255 
256 inline lsst::geom::Point2D computeSrcPos(
257  int destCol,
258  int destRow,
259  lsst::geom::Point2D const &destXY0,
260  geom::SkyWcs const &destWcs,
261  geom::SkyWcs const &srcWcs)
262 {
263  lsst::geom::Point2D const destPix(image::indexToPosition(destCol + destXY0[0]),
264  image::indexToPosition(destRow + destXY0[1]));
265  return srcWcs.skyToPixel(destWcs.pixelToSky(destPix));
266 }
267 
268 inline double computeRelativeArea(
269  lsst::geom::Point2D const &srcPos,
270  lsst::geom::Point2D const
271  &leftSrcPos,
272  lsst::geom::Point2D const &upSrcPos)
273 {
274  lsst::geom::Extent2D dSrcA = srcPos - leftSrcPos;
275  lsst::geom::Extent2D dSrcB = srcPos - upSrcPos;
276 
277  return std::abs(dSrcA.getX() * dSrcB.getY() - dSrcA.getY() * dSrcB.getX());
278 }
279 
280 } // namespace
281 
282 template <typename DestImageT, typename SrcImageT>
283 int warpImage(DestImageT &destImage, geom::SkyWcs const &destWcs, SrcImageT const &srcImage,
284  geom::SkyWcs const &srcWcs, WarpingControl const &control,
285  typename DestImageT::SinglePixel padValue) {
286  auto srcToDest = geom::makeWcsPairTransform(srcWcs, destWcs);
287  return warpImage(destImage, srcImage, *srcToDest, control, padValue);
288 }
289 
290 template <typename DestImageT, typename SrcImageT>
291 int warpImage(DestImageT &destImage, SrcImageT const &srcImage,
292  geom::TransformPoint2ToPoint2 const &srcToDest, WarpingControl const &control,
293  typename DestImageT::SinglePixel padValue) {
294  if (imagesOverlap(destImage, srcImage)) {
295  throw LSST_EXCEPT(pexExcept::InvalidParameterError, "destImage overlaps srcImage; cannot warp");
296  }
297  if (destImage.getBBox(image::LOCAL).isEmpty()) {
298  return 0;
299  }
300  // if src image is too small then don't try to warp
301  std::shared_ptr<SeparableKernel> warpingKernelPtr = control.getWarpingKernel();
302  try {
303  warpingKernelPtr->shrinkBBox(srcImage.getBBox(image::LOCAL));
305  for (int y = 0, height = destImage.getHeight(); y < height; ++y) {
306  for (typename DestImageT::x_iterator destPtr = destImage.row_begin(y), end = destImage.row_end(y);
307  destPtr != end; ++destPtr) {
308  *destPtr = padValue;
309  }
310  }
311  return 0;
312  }
313  int interpLength = control.getInterpLength();
314 
315  std::shared_ptr<LanczosWarpingKernel const> const lanczosKernelPtr =
316  std::dynamic_pointer_cast<LanczosWarpingKernel>(warpingKernelPtr);
317 
318  int numGoodPixels = 0;
319 
320  // compute a transform from local destination pixels to parent source pixels
321  auto const parentDestToParentSrc = srcToDest.inverted();
322  std::vector<double> const localDestToParentDestVec = {static_cast<double>(destImage.getX0()),
323  static_cast<double>(destImage.getY0())};
324  auto const localDestToParentDest = geom::TransformPoint2ToPoint2(ast::ShiftMap(localDestToParentDestVec));
325  auto const localDestToParentSrc = localDestToParentDest.then(*parentDestToParentSrc);
326 
327  // Get the source MaskedImage and a pixel accessor to it.
328  int const srcWidth = srcImage.getWidth();
329  int const srcHeight = srcImage.getHeight();
330  LOGL_DEBUG("TRACE2.afw.math.warp", "source image width=%d; height=%d", srcWidth, srcHeight);
331 
332  int const destWidth = destImage.getWidth();
333  int const destHeight = destImage.getHeight();
334  LOGL_DEBUG("TRACE2.afw.math.warp", "remap image width=%d; height=%d", destWidth, destHeight);
335 
336  // Set each pixel of destExposure's MaskedImage
337  LOGL_DEBUG("TRACE3.afw.math.warp", "Remapping masked image");
338 
339  int const maxCol = destWidth - 1;
340  int const maxRow = destHeight - 1;
341 
342  detail::WarpAtOnePoint<DestImageT, SrcImageT> warpAtOnePoint(srcImage, control, padValue);
343 
344  if (interpLength > 0) {
345  // Use interpolation. Note that 1 produces the same result as no interpolation
346  // but uses this code branch, thus providing an easy way to compare the two branches.
347 
348  // Estimate for number of horizontal interpolation band edges, to reserve memory in vectors
349  int const numColEdges = 2 + ((destWidth - 1) / interpLength);
350 
351  // A list of edge column indices for interpolation bands;
352  // starts at -1, increments by interpLen (except the final interval), and ends at destWidth-1
353  std::vector<int> edgeColList;
354  edgeColList.reserve(numColEdges);
355 
356  // A list of 1/column width for horizontal interpolation bands; the first value is garbage.
357  // The inverse is used for speed because the values are always multiplied.
358  std::vector<double> invWidthList;
359  invWidthList.reserve(numColEdges);
360 
361  // Compute edgeColList and invWidthList
362  edgeColList.push_back(-1);
363  invWidthList.push_back(0.0);
364  for (int prevEndCol = -1; prevEndCol < maxCol; prevEndCol += interpLength) {
365  int endCol = prevEndCol + interpLength;
366  if (endCol > maxCol) {
367  endCol = maxCol;
368  }
369  edgeColList.push_back(endCol);
370  assert(endCol - prevEndCol > 0);
371  invWidthList.push_back(1.0 / static_cast<double>(endCol - prevEndCol));
372  }
373  assert(edgeColList.back() == maxCol);
374 
375  // A list of delta source positions along the edge columns of the horizontal interpolation bands
376  std::vector<lsst::geom::Extent2D> yDeltaSrcPosList(edgeColList.size());
377 
378  // A cache of pixel positions on the source corresponding to the previous or current row
379  // of the destination image.
380  // The first value is for column -1 because the previous source position is used to compute relative
381  // area To simplify the indexing, use an iterator that starts at begin+1, thus: srcPosView =
382  // srcPosList.begin() + 1 srcPosView[col-1] and lower indices are for this row srcPosView[col] and
383  // higher indices are for the previous row
384  std::vector<lsst::geom::Point2D> srcPosList(1 + destWidth);
385  std::vector<lsst::geom::Point2D>::iterator const srcPosView = srcPosList.begin() + 1;
386 
387  std::vector<lsst::geom::Point2D> endColPosList;
388  endColPosList.reserve(numColEdges);
389 
390  // Initialize srcPosList for row -1
391  for (int colBand = 0, endBand = edgeColList.size(); colBand < endBand; ++colBand) {
392  int const endCol = edgeColList[colBand];
393  endColPosList.emplace_back(lsst::geom::Point2D(endCol, -1));
394  }
395  auto rightSrcPosList = localDestToParentSrc->applyForward(endColPosList);
396  srcPosView[-1] = rightSrcPosList[0];
397  for (int colBand = 1, endBand = edgeColList.size(); colBand < endBand; ++colBand) {
398  int const prevEndCol = edgeColList[colBand - 1];
399  int const endCol = edgeColList[colBand];
400  lsst::geom::Point2D leftSrcPos = srcPosView[prevEndCol];
401 
402  lsst::geom::Extent2D xDeltaSrcPos =
403  (rightSrcPosList[colBand] - leftSrcPos) * invWidthList[colBand];
404 
405  for (int col = prevEndCol + 1; col <= endCol; ++col) {
406  srcPosView[col] = srcPosView[col - 1] + xDeltaSrcPos;
407  }
408  }
409 
410  int endRow = -1;
411  while (endRow < maxRow) {
412  // Next horizontal interpolation band
413 
414  int prevEndRow = endRow;
415  endRow = prevEndRow + interpLength;
416  if (endRow > maxRow) {
417  endRow = maxRow;
418  }
419  assert(endRow - prevEndRow > 0);
420  double interpInvHeight = 1.0 / static_cast<double>(endRow - prevEndRow);
421 
422  // Set yDeltaSrcPosList for this horizontal interpolation band
423  std::vector<lsst::geom::Point2D> destRowPosList;
424  destRowPosList.reserve(edgeColList.size());
425  for (int colBand = 0, endBand = edgeColList.size(); colBand < endBand; ++colBand) {
426  int endCol = edgeColList[colBand];
427  destRowPosList.emplace_back(lsst::geom::Point2D(endCol, endRow));
428  }
429  auto bottomSrcPosList = localDestToParentSrc->applyForward(destRowPosList);
430  for (int colBand = 0, endBand = edgeColList.size(); colBand < endBand; ++colBand) {
431  int endCol = edgeColList[colBand];
432  yDeltaSrcPosList[colBand] =
433  (bottomSrcPosList[colBand] - srcPosView[endCol]) * interpInvHeight;
434  }
435 
436  for (int row = prevEndRow + 1; row <= endRow; ++row) {
437  typename DestImageT::x_iterator destXIter = destImage.row_begin(row);
438  srcPosView[-1] += yDeltaSrcPosList[0];
439  for (int colBand = 1, endBand = edgeColList.size(); colBand < endBand; ++colBand) {
440  // Next vertical interpolation band
441 
442  int const prevEndCol = edgeColList[colBand - 1];
443  int const endCol = edgeColList[colBand];
444 
445  // Compute xDeltaSrcPos; remember that srcPosView contains
446  // positions for this row in prevEndCol and smaller indices,
447  // and positions for the previous row for larger indices (including endCol)
448  lsst::geom::Point2D leftSrcPos = srcPosView[prevEndCol];
449  lsst::geom::Point2D rightSrcPos = srcPosView[endCol] + yDeltaSrcPosList[colBand];
450  lsst::geom::Extent2D xDeltaSrcPos = (rightSrcPos - leftSrcPos) * invWidthList[colBand];
451 
452  for (int col = prevEndCol + 1; col <= endCol; ++col, ++destXIter) {
453  lsst::geom::Point2D leftSrcPos = srcPosView[col - 1];
454  lsst::geom::Point2D srcPos = leftSrcPos + xDeltaSrcPos;
455  double relativeArea = computeRelativeArea(srcPos, leftSrcPos, srcPosView[col]);
456 
457  srcPosView[col] = srcPos;
458 
459  if (warpAtOnePoint(
460  destXIter, srcPos, relativeArea,
462  ++numGoodPixels;
463  }
464  } // for col
465  } // for col band
466  } // for row
467  } // while next row band
468 
469  } else {
470  // No interpolation
471 
472  // prevSrcPosList = source positions from the previous row; these are used to compute pixel area;
473  // to begin, compute sources positions corresponding to destination row = -1
475  destPosList.reserve(1 + destWidth);
476  for (int col = -1; col < destWidth; ++col) {
477  destPosList.emplace_back(lsst::geom::Point2D(col, -1));
478  }
479  auto prevSrcPosList = localDestToParentSrc->applyForward(destPosList);
480 
481  for (int row = 0; row < destHeight; ++row) {
482  destPosList.clear();
483  for (int col = -1; col < destWidth; ++col) {
484  destPosList.emplace_back(lsst::geom::Point2D(col, row));
485  }
486  auto srcPosList = localDestToParentSrc->applyForward(destPosList);
487 
488  typename DestImageT::x_iterator destXIter = destImage.row_begin(row);
489  for (int col = 0; col < destWidth; ++col, ++destXIter) {
490  // column index = column + 1 because the first entry in srcPosList is for column -1
491  auto srcPos = srcPosList[col + 1];
492  double relativeArea =
493  computeRelativeArea(srcPos, prevSrcPosList[col], prevSrcPosList[col + 1]);
494 
495  if (warpAtOnePoint(destXIter, srcPos, relativeArea,
497  ++numGoodPixels;
498  }
499  } // for col
500  // move points from srcPosList to prevSrcPosList (we don't care about what ends up in srcPosList
501  // because it will be reallocated anyway)
502  swap(srcPosList, prevSrcPosList);
503  } // for row
504  } // if interp
505 
506  return numGoodPixels;
507 }
508 
509 template <typename DestImageT, typename SrcImageT>
510 int warpCenteredImage(DestImageT &destImage, SrcImageT const &srcImage,
511  lsst::geom::LinearTransform const &linearTransform,
512  lsst::geom::Point2D const &centerPosition, WarpingControl const &control,
513  typename DestImageT::SinglePixel padValue) {
514  // force src and dest to be the same size and xy0
515  if ((destImage.getWidth() != srcImage.getWidth()) || (destImage.getHeight() != srcImage.getHeight()) ||
516  (destImage.getXY0() != srcImage.getXY0())) {
517  std::ostringstream errStream;
518  errStream << "src and dest images must have same size and xy0.";
520  }
521 
522  // set the xy0 coords to 0,0 to make life easier
523  SrcImageT srcImageCopy(srcImage, true);
524  srcImageCopy.setXY0(0, 0);
525  destImage.setXY0(0, 0);
526  lsst::geom::Extent2D cLocal =
527  lsst::geom::Extent2D(centerPosition) - lsst::geom::Extent2D(srcImage.getXY0());
528 
529  // for the affine transform, the centerPosition will not only get sheared, but also
530  // moved slightly. So we'll include a translation to move it back by an amount
531  // centerPosition - translatedCenterPosition
532  lsst::geom::AffineTransform affTran(linearTransform, cLocal - linearTransform(cLocal));
534 
535 // now warp
536 #if 0
537  static float t = 0.0;
538  float t_before = 1.0*clock()/CLOCKS_PER_SEC;
539  int n = warpImage(destImage, srcImageCopy, affTran, control, padValue);
540  float t_after = 1.0*clock()/CLOCKS_PER_SEC;
541  float dt = t_after - t_before;
542  t += dt;
543  std::cout <<srcImage.getWidth()<<"x"<<srcImage.getHeight()<<": "<< dt <<" "<< t <<std::endl;
544 #else
545  int n = warpImage(destImage, srcImageCopy, *affineTransform22, control, padValue);
546 #endif
547 
548  // fix the origin and we're done.
549  destImage.setXY0(srcImage.getXY0());
550 
551  return n;
552 }
553 
554 //
555 // Explicit instantiations
556 //
558 // may need to omit default params for EXPOSURE -- original code did that and it worked
559 #define EXPOSURE(PIXTYPE) image::Exposure<PIXTYPE, image::MaskPixel, image::VariancePixel>
560 #define MASKEDIMAGE(PIXTYPE) image::MaskedImage<PIXTYPE, image::MaskPixel, image::VariancePixel>
561 #define IMAGE(PIXTYPE) image::Image<PIXTYPE>
562 #define NL /* */
563 
564 #define INSTANTIATE(DESTIMAGEPIXELT, SRCIMAGEPIXELT) \
565  template int warpCenteredImage( \
566  IMAGE(DESTIMAGEPIXELT) & destImage, IMAGE(SRCIMAGEPIXELT) const &srcImage, \
567  lsst::geom::LinearTransform const &linearTransform, lsst::geom::Point2D const &centerPosition, \
568  WarpingControl const &control, IMAGE(DESTIMAGEPIXELT)::SinglePixel padValue); \
569  NL template int warpCenteredImage( \
570  MASKEDIMAGE(DESTIMAGEPIXELT) & destImage, MASKEDIMAGE(SRCIMAGEPIXELT) const &srcImage, \
571  lsst::geom::LinearTransform const &linearTransform, lsst::geom::Point2D const &centerPosition, \
572  WarpingControl const &control, MASKEDIMAGE(DESTIMAGEPIXELT)::SinglePixel padValue); \
573  NL template int warpImage(IMAGE(DESTIMAGEPIXELT) & destImage, IMAGE(SRCIMAGEPIXELT) const &srcImage, \
574  geom::TransformPoint2ToPoint2 const &srcToDest, WarpingControl const &control, \
575  IMAGE(DESTIMAGEPIXELT)::SinglePixel padValue); \
576  NL template int warpImage(MASKEDIMAGE(DESTIMAGEPIXELT) & destImage, \
577  MASKEDIMAGE(SRCIMAGEPIXELT) const &srcImage, \
578  geom::TransformPoint2ToPoint2 const &srcToDest, WarpingControl const &control, \
579  MASKEDIMAGE(DESTIMAGEPIXELT)::SinglePixel padValue); \
580  NL template int warpImage(IMAGE(DESTIMAGEPIXELT) & destImage, geom::SkyWcs const &destWcs, \
581  IMAGE(SRCIMAGEPIXELT) const &srcImage, geom::SkyWcs const &srcWcs, \
582  WarpingControl const &control, IMAGE(DESTIMAGEPIXELT)::SinglePixel padValue); \
583  NL template int warpImage(MASKEDIMAGE(DESTIMAGEPIXELT) & destImage, geom::SkyWcs const &destWcs, \
584  MASKEDIMAGE(SRCIMAGEPIXELT) const &srcImage, geom::SkyWcs const &srcWcs, \
585  WarpingControl const &control, \
586  MASKEDIMAGE(DESTIMAGEPIXELT)::SinglePixel padValue); \
587  NL template int warpExposure(EXPOSURE(DESTIMAGEPIXELT) & destExposure, \
588  EXPOSURE(SRCIMAGEPIXELT) const &srcExposure, WarpingControl const &control, \
589  EXPOSURE(DESTIMAGEPIXELT)::MaskedImageT::SinglePixel padValue);
590 
591 INSTANTIATE(double, double)
592 INSTANTIATE(double, float)
593 INSTANTIATE(double, int)
594 INSTANTIATE(double, std::uint16_t)
595 INSTANTIATE(float, float)
596 INSTANTIATE(float, int)
598 INSTANTIATE(int, int)
601 } // namespace math
602 } // namespace afw
603 } // namespace lsst
table::Key< std::string > name
Definition: Amplifier.cc:116
int end
double x
#define INSTANTIATE(FROMSYS, TOSYS)
Definition: Detector.cc:484
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
LSST DM logging module built on log4cxx.
#define LOGL_DEBUG(logger, message...)
Log a debug-level message using a varargs/printf style interface.
Definition: Log.h:504
Implementation of the Photometric Calibration class.
std::ostream * os
Definition: Schema.cc:746
std::string prefix
Definition: SchemaMapper.cc:79
int y
Definition: SpanSet.cc:49
T back(T... args)
T begin(T... args)
ShiftMap is a linear Mapping which shifts each axis by a specified constant value.
Definition: ShiftMap.h:40
std::string toString(std::string const &="") const override
Return string representation.
Kernel::Pixel operator()(double x) const override
Solve bilinear equation.
std::shared_ptr< Kernel > clone() const override
Return a pointer to a deep copy of this kernel.
void setKernelParameter(unsigned int ind, double value) const override
Set one kernel parameter.
std::vector< double > _params
Definition: Function.h:185
lsst::geom::Extent2I const getDimensions() const
Return the Kernel's dimensions (width, height)
Definition: Kernel.h:213
lsst::geom::Point2I getCtr() const
Return index of kernel's center.
Definition: Kernel.h:235
int getWidth() const
Return the Kernel's width.
Definition: Kernel.h:225
Lanczos warping: accurate but slow and can introduce ringing artifacts.
Definition: warpExposure.h:65
int getOrder() const
get the order of the kernel
void setKernelParameter(unsigned int ind, double value) const override
Set one kernel parameter.
std::shared_ptr< Kernel > clone() const override
Return a pointer to a deep copy of this kernel.
Definition: warpExposure.cc:97
Kernel::Pixel operator()(double x) const override
Solve nearest neighbor equation.
std::string toString(std::string const &="") const override
Return string representation.
Nearest neighbor warping: fast; good for undersampled data.
Definition: warpExposure.h:163
std::shared_ptr< Kernel > clone() const override
Return a pointer to a deep copy of this kernel.
void setKernelParameter(unsigned int ind, double value) const override
Set one kernel parameter.
A kernel described by a pair of functions: func(x, y) = colFunc(x) * rowFunc(y)
Definition: Kernel.h:861
std::shared_ptr< Kernel > clone() const override
Return a pointer to a deep copy of this kernel.
void setKernelParameter(unsigned int ind, double value) const override
Set one kernel parameter.
Parameters to control convolution.
Definition: warpExposure.h:250
void setWarpingKernel(SeparableKernel const &warpingKernel)
set the warping kernel
int getInterpLength() const
get the interpolation length (pixels)
Definition: warpExposure.h:303
void setWarpingKernelName(std::string const &warpingKernelName)
set the warping kernel by name
void setMaskWarpingKernelName(std::string const &maskWarpingKernelName)
set or clear the mask warping kernel by name
void setMaskWarpingKernel(SeparableKernel const &maskWarpingKernel)
set the mask warping kernel
std::shared_ptr< SeparableKernel > getWarpingKernel() const
get the warping kernel
std::shared_ptr< SeparableKernel > getMaskWarpingKernel() const
get the mask warping kernel
A functor that computes one warped pixel.
An affine coordinate transformation consisting of a linear transformation and an offset.
An integer coordinate rectangle.
Definition: Box.h:55
bool contains(Point2I const &point) const noexcept
Return true if the box contains the point.
Definition: Box.cc:114
A 2D linear coordinate transformation.
Reports invalid arguments.
Definition: Runtime.h:66
T clear(T... args)
T emplace_back(T... args)
T empty(T... args)
T endl(T... args)
void swap(CameraSys &a, CameraSys &b)
Definition: CameraSys.h:157
std::shared_ptr< TransformPoint2ToPoint2 > makeTransform(lsst::geom::AffineTransform const &affine)
Wrap an lsst::geom::AffineTransform as a Transform.
std::shared_ptr< TransformPoint2ToPoint2 > makeWcsPairTransform(SkyWcs const &src, SkyWcs const &dst)
A Transform obtained by putting two SkyWcs objects "back to back".
Definition: SkyWcs.cc:151
Transform< Point2Endpoint, Point2Endpoint > TransformPoint2ToPoint2
Definition: Transform.h:300
double indexToPosition(double ind)
Convert image index to image position.
Definition: ImageUtils.h:55
bool imagesOverlap(ImageBase< T1 > const &image1, ImageBase< T2 > const &image2)
Return true if the pixels for two images or masks overlap in memory.
Definition: Image.cc:706
std::shared_ptr< SeparableKernel > makeWarpingKernel(std::string name)
Return a warping kernel given its name.
int warpCenteredImage(DestImageT &destImage, SrcImageT const &srcImage, lsst::geom::LinearTransform const &linearTransform, lsst::geom::Point2D const &centerPosition, WarpingControl const &control, typename DestImageT::SinglePixel padValue=lsst::afw::math::edgePixel< DestImageT >(typename lsst::afw::image::detail::image_traits< DestImageT >::image_category()))
Warp an image with a LinearTranform about a specified point.
int warpImage(DestImageT &destImage, geom::SkyWcs const &destWcs, SrcImageT const &srcImage, geom::SkyWcs const &srcWcs, WarpingControl const &control, typename DestImageT::SinglePixel padValue=lsst::afw::math::edgePixel< DestImageT >(typename lsst::afw::image::detail::image_traits< DestImageT >::image_category()))
Warp an Image or MaskedImage to a new Wcs.
int warpExposure(DestExposureT &destExposure, SrcExposureT const &srcExposure, WarpingControl const &control, typename DestExposureT::MaskedImageT::SinglePixel padValue=lsst::afw::math::edgePixel< typename DestExposureT::MaskedImageT >(typename lsst::afw::image::detail::image_traits< typename DestExposureT::MaskedImageT >::image_category()))
Warp (remap) one exposure to another.
Extent< double, 2 > Extent2D
Definition: Extent.h:400
Angle abs(Angle const &a)
Definition: Angle.h:106
A base class for image defects.
T push_back(T... args)
T reserve(T... args)
T size(T... args)
afw::table::Key< std::string > warpingKernelName
Definition: CoaddPsf.cc:341
int row
Definition: CR.cc:145
int col
Definition: CR.cc:144
T str(T... args)
ImageT::image_category image_category
Definition: ImageBase.h:67
T swap(T... args)