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
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 <limits>
31 #include <memory>
32 #include <sstream>
33 #include <string>
34 #include <vector>
35 #include <utility>
36 #include <regex>
37 
38 #include "astshim.h"
39 
40 #include "lsst/log/Log.h"
41 #include "lsst/pex/exceptions.h"
42 #include "lsst/geom.h"
44 #include "lsst/afw/geom.h"
45 #include "lsst/afw/math/Kernel.h"
51 #include "lsst/afw/table/io/Persistable.cc" // Needed for PersistableFacade::dynamicCast
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 
161 namespace {
162 
163 struct LanczosKernelPersistenceHelper {
164  table::Schema schema;
165  table::Key<int> order;
166 
167  static LanczosKernelPersistenceHelper const &get() {
168  static LanczosKernelPersistenceHelper const instance;
169  return instance;
170  }
171 
172  LanczosKernelPersistenceHelper(LanczosKernelPersistenceHelper const &) = delete;
173  LanczosKernelPersistenceHelper(LanczosKernelPersistenceHelper &&) = delete;
174  LanczosKernelPersistenceHelper &operator=(LanczosKernelPersistenceHelper const &) = delete;
175  LanczosKernelPersistenceHelper &operator=(LanczosKernelPersistenceHelper &&) = delete;
176 
177 private:
178  LanczosKernelPersistenceHelper()
179  : schema(), order(schema.addField<int>("order", "order of Lanczos function")) {}
180 };
181 
182 class : public table::io::PersistableFactory {
183  std::shared_ptr<table::io::Persistable> read(table::io::InputArchive const &archive,
184  table::io::CatalogVector const &catalogs) const override {
185  auto const &keys = LanczosKernelPersistenceHelper::get();
186  LSST_ARCHIVE_ASSERT(catalogs.size() == 1u);
187  LSST_ARCHIVE_ASSERT(catalogs.front().size() == 1u);
188  afw::table::BaseRecord const &record = catalogs.front().front();
189  LSST_ARCHIVE_ASSERT(record.getSchema() == keys.schema);
190  return std::make_shared<LanczosWarpingKernel>(record.get(keys.order));
191  }
192 
194 } lanczosFactory("LanczosWarpingKernel");
195 
196 template <class T>
197 class DefaultPersistableFactory : public table::io::PersistableFactory {
198  std::shared_ptr<table::io::Persistable> read(table::io::InputArchive const &archive,
199  table::io::CatalogVector const &catalogs) const override {
200  LSST_ARCHIVE_ASSERT(catalogs.empty());
201  return std::make_shared<T>();
202  }
203 
205 };
206 
207 DefaultPersistableFactory<BilinearWarpingKernel> bilinearFactory("BilinearWarpingKernel");
208 DefaultPersistableFactory<NearestWarpingKernel> nearestFactory("NearestWarpingKernel");
209 
210 } // namespace
211 
213  auto const &keys = LanczosKernelPersistenceHelper::get();
214  table::BaseCatalog catalog = handle.makeCatalog(keys.schema);
215  std::shared_ptr<table::BaseRecord> record = catalog.addNew();
216  record->set(keys.order, getOrder());
217  handle.saveCatalog(catalog);
218 }
219 
221 
223 
225  using KernelPtr = std::shared_ptr<SeparableKernel>;
226  std::smatch matches;
227  static const std::regex LanczosRE("lanczos(\\d+)");
228  if (name == "bilinear") {
229  return KernelPtr(new BilinearWarpingKernel());
230  } else if (std::regex_match(name, matches, LanczosRE)) {
231  std::string orderStr(matches[1].first, matches[1].second);
232  int order = std::stoi(orderStr);
233  return KernelPtr(new LanczosWarpingKernel(order));
234  } else if (name == "nearest") {
235  return KernelPtr(new NearestWarpingKernel());
236  } else {
237  throw LSST_EXCEPT(pexExcept::InvalidParameterError, "unknown warping kernel name: \"" + name + "\"");
238  }
239 }
240 
242  if (_warpingKernelPtr->getCacheSize() != _cacheSize) {
243  _warpingKernelPtr->computeCache(_cacheSize);
244  }
245  return _warpingKernelPtr;
246 };
247 
250  setWarpingKernel(*warpingKernelPtr);
251 }
252 
254  if (_maskWarpingKernelPtr) {
255  _testWarpingKernels(warpingKernel, *_maskWarpingKernelPtr);
256  }
257  std::shared_ptr<SeparableKernel> warpingKernelPtr(
258  std::static_pointer_cast<SeparableKernel>(warpingKernel.clone()));
259  _warpingKernelPtr = warpingKernelPtr;
260 }
261 
263  if (_maskWarpingKernelPtr) { // lazily update kernel cache
264  if (_maskWarpingKernelPtr->getCacheSize() != _cacheSize) {
265  _maskWarpingKernelPtr->computeCache(_cacheSize);
266  }
267  }
268  return _maskWarpingKernelPtr;
269 }
270 
271 void WarpingControl::setMaskWarpingKernelName(std::string const &maskWarpingKernelName) {
272  if (!maskWarpingKernelName.empty()) {
273  std::shared_ptr<SeparableKernel> maskWarpingKernelPtr(makeWarpingKernel(maskWarpingKernelName));
274  setMaskWarpingKernel(*maskWarpingKernelPtr);
275  } else {
276  _maskWarpingKernelPtr.reset();
277  }
278 }
279 
281  _testWarpingKernels(*_warpingKernelPtr, maskWarpingKernel);
282  _maskWarpingKernelPtr = std::static_pointer_cast<SeparableKernel>(maskWarpingKernel.clone());
283 }
284 
285 void WarpingControl::_testWarpingKernels(SeparableKernel const &warpingKernel,
286  SeparableKernel const &maskWarpingKernel) const {
287  lsst::geom::Box2I kernelBBox =
289  warpingKernel.getDimensions());
290  lsst::geom::Box2I maskKernelBBox =
292  maskWarpingKernel.getDimensions());
293  if (!kernelBBox.contains(maskKernelBBox)) {
295  "warping kernel is smaller than mask warping kernel");
296  }
297 }
298 
299 namespace {
300 
301 struct WarpingControlPersistenceHelper {
303  table::Key<int> warpingKernelIndex;
304  table::Key<table::Flag> hasMaskKernel;
305  table::Key<int> maskKernelIndex; // undefined if !hasMaskKernel
306  table::Key<int> cacheSize;
307  table::Key<int> interpLength;
308  table::Key<image::MaskPixel> growFullMask;
309 
310  static WarpingControlPersistenceHelper const &get() {
311  static WarpingControlPersistenceHelper const instance;
312  return instance;
313  }
314 
315  WarpingControlPersistenceHelper(WarpingControlPersistenceHelper const &) = delete;
316  WarpingControlPersistenceHelper(WarpingControlPersistenceHelper &&) = delete;
317  WarpingControlPersistenceHelper &operator=(WarpingControlPersistenceHelper const &) = delete;
318  WarpingControlPersistenceHelper &operator=(WarpingControlPersistenceHelper &&) = delete;
319 
320 private:
321  WarpingControlPersistenceHelper()
322  : schema(),
324  schema.addField<int>("warpingKernelIndex", "archive ID of nested warping kernel")),
325  hasMaskKernel(schema.addField<table::Flag>("hasMaskKernel", "whether a mask kernel is stored")),
326  maskKernelIndex(schema.addField<int>("maskKernelIndex",
327  "archive ID of nested mask kernel. "
328  "Valid only if hasMaskKernel")),
329  cacheSize(schema.addField<int>("cacheSize", "Cache size for warping kernel(s)")),
330  interpLength(schema.addField<int>("interpLength",
331  "Distance over which WCS can be linearly interpolated")),
332  growFullMask(schema.addField<image::MaskPixel>(
333  "growFullMask", "bits to grow to full width of image/variance kernel")) {}
334 };
335 
336 std::string _getWarpingControlPersistenceName() { return "WarpingControl"; }
337 
338 class : public table::io::PersistableFactory {
339  std::shared_ptr<table::io::Persistable> read(table::io::InputArchive const &archive,
340  table::io::CatalogVector const &catalogs) const override {
341  auto const &keys = WarpingControlPersistenceHelper::get();
342  LSST_ARCHIVE_ASSERT(catalogs.size() == 1u);
343  LSST_ARCHIVE_ASSERT(catalogs.front().size() == 1u);
344  afw::table::BaseRecord const &record = catalogs.front().front();
345  LSST_ARCHIVE_ASSERT(record.getSchema() == keys.schema);
346 
347  // Warping kernels are dummy values, set true kernels later
348  auto control = std::make_shared<WarpingControl>("bilinear", "", record.get(keys.cacheSize),
349  record.get(keys.interpLength),
350  record.get(keys.growFullMask));
351  // archive.get returns a shared_ptr, so this code depends on the
352  // undocumented fact that setWarpingKernel and setMaskWarpingKernel
353  // make defensive copies.
354  control->setWarpingKernel(*archive.get<SeparableKernel>(record.get(keys.warpingKernelIndex)));
355  if (record.get(keys.hasMaskKernel)) {
356  control->setMaskWarpingKernel(*archive.get<SeparableKernel>(record.get(keys.maskKernelIndex)));
357  }
358  return control;
359  }
360 
362 } controlFactory(_getWarpingControlPersistenceName());
363 
364 } // namespace
365 
366 std::string WarpingControl::getPersistenceName() const { return _getWarpingControlPersistenceName(); }
367 
368 std::string WarpingControl::getPythonModule() const { return "lsst.afw.math"; }
369 
370 bool WarpingControl::isPersistable() const noexcept {
371  return _warpingKernelPtr->isPersistable() &&
372  (!hasMaskWarpingKernel() || _maskWarpingKernelPtr->isPersistable());
373 }
374 
376  auto const &keys = WarpingControlPersistenceHelper::get();
377  table::BaseCatalog catalog = handle.makeCatalog(keys.schema);
378  std::shared_ptr<table::BaseRecord> record = catalog.addNew();
379 
380  record->set(keys.warpingKernelIndex, handle.put(_warpingKernelPtr));
381  record->set(keys.hasMaskKernel, hasMaskWarpingKernel());
382  if (hasMaskWarpingKernel()) {
383  record->set(keys.maskKernelIndex, handle.put(_maskWarpingKernelPtr));
384  }
385  record->set(keys.cacheSize, _cacheSize);
386  record->set(keys.interpLength, _interpLength);
387  record->set(keys.growFullMask, _growFullMask);
388 
389  handle.saveCatalog(catalog);
390 }
391 
392 template <typename DestExposureT, typename SrcExposureT>
393 int warpExposure(DestExposureT &destExposure, SrcExposureT const &srcExposure, WarpingControl const &control,
394  typename DestExposureT::MaskedImageT::SinglePixel padValue) {
395  if (!destExposure.hasWcs()) {
396  throw LSST_EXCEPT(pexExcept::InvalidParameterError, "destExposure has no Wcs");
397  }
398  if (!srcExposure.hasWcs()) {
399  throw LSST_EXCEPT(pexExcept::InvalidParameterError, "srcExposure has no Wcs");
400  }
401  typename DestExposureT::MaskedImageT mi = destExposure.getMaskedImage();
402  if (srcExposure.getInfo()->hasId()) {
403  destExposure.getInfo()->setId(srcExposure.getInfo()->getId());
404  }
405  destExposure.setPhotoCalib(srcExposure.getPhotoCalib());
406  destExposure.setFilterLabel(srcExposure.getFilterLabel());
407  destExposure.getInfo()->setVisitInfo(srcExposure.getInfo()->getVisitInfo());
408  return warpImage(mi, *destExposure.getWcs(), srcExposure.getMaskedImage(), *srcExposure.getWcs(), control,
409  padValue);
410 }
411 
412 namespace {
413 
414 inline lsst::geom::Point2D computeSrcPos(
415  int destCol,
416  int destRow,
417  lsst::geom::Point2D const &destXY0,
418  geom::SkyWcs const &destWcs,
419  geom::SkyWcs const &srcWcs)
420 {
421  lsst::geom::Point2D const destPix(image::indexToPosition(destCol + destXY0[0]),
422  image::indexToPosition(destRow + destXY0[1]));
423  return srcWcs.skyToPixel(destWcs.pixelToSky(destPix));
424 }
425 
426 inline double computeRelativeArea(
427  lsst::geom::Point2D const &srcPos,
428  lsst::geom::Point2D const
429  &leftSrcPos,
430  lsst::geom::Point2D const &upSrcPos)
431 {
432  lsst::geom::Extent2D dSrcA = srcPos - leftSrcPos;
433  lsst::geom::Extent2D dSrcB = srcPos - upSrcPos;
434 
435  return std::abs(dSrcA.getX() * dSrcB.getY() - dSrcA.getY() * dSrcB.getX());
436 }
437 
438 } // namespace
439 
440 template <typename DestImageT, typename SrcImageT>
441 int warpImage(DestImageT &destImage, geom::SkyWcs const &destWcs, SrcImageT const &srcImage,
442  geom::SkyWcs const &srcWcs, WarpingControl const &control,
443  typename DestImageT::SinglePixel padValue) {
444  auto srcToDest = geom::makeWcsPairTransform(srcWcs, destWcs);
445  return warpImage(destImage, srcImage, *srcToDest, control, padValue);
446 }
447 
448 template <typename DestImageT, typename SrcImageT>
449 int warpImage(DestImageT &destImage, SrcImageT const &srcImage,
450  geom::TransformPoint2ToPoint2 const &srcToDest, WarpingControl const &control,
451  typename DestImageT::SinglePixel padValue) {
452  if (imagesOverlap(destImage, srcImage)) {
453  throw LSST_EXCEPT(pexExcept::InvalidParameterError, "destImage overlaps srcImage; cannot warp");
454  }
455  if (destImage.getBBox(image::LOCAL).isEmpty()) {
456  return 0;
457  }
458  // if src image is too small then don't try to warp
459  std::shared_ptr<SeparableKernel> warpingKernelPtr = control.getWarpingKernel();
460  try {
461  warpingKernelPtr->shrinkBBox(srcImage.getBBox(image::LOCAL));
463  for (int y = 0, height = destImage.getHeight(); y < height; ++y) {
464  for (typename DestImageT::x_iterator destPtr = destImage.row_begin(y), end = destImage.row_end(y);
465  destPtr != end; ++destPtr) {
466  *destPtr = padValue;
467  }
468  }
469  return 0;
470  }
471  int interpLength = control.getInterpLength();
472 
473  std::shared_ptr<LanczosWarpingKernel const> const lanczosKernelPtr =
474  std::dynamic_pointer_cast<LanczosWarpingKernel>(warpingKernelPtr);
475 
476  int numGoodPixels = 0;
477 
478  // compute a transform from local destination pixels to parent source pixels
479  auto const parentDestToParentSrc = srcToDest.inverted();
480  std::vector<double> const localDestToParentDestVec = {static_cast<double>(destImage.getX0()),
481  static_cast<double>(destImage.getY0())};
482  auto const localDestToParentDest = geom::TransformPoint2ToPoint2(ast::ShiftMap(localDestToParentDestVec));
483  auto const localDestToParentSrc = localDestToParentDest.then(*parentDestToParentSrc);
484 
485  // Get the source MaskedImage and a pixel accessor to it.
486  int const srcWidth = srcImage.getWidth();
487  int const srcHeight = srcImage.getHeight();
488  LOGL_DEBUG("TRACE2.lsst.afw.math.warp", "source image width=%d; height=%d", srcWidth, srcHeight);
489 
490  int const destWidth = destImage.getWidth();
491  int const destHeight = destImage.getHeight();
492  LOGL_DEBUG("TRACE2.lsst.afw.math.warp", "remap image width=%d; height=%d", destWidth, destHeight);
493 
494  // Set each pixel of destExposure's MaskedImage
495  LOGL_DEBUG("TRACE3.lsst.afw.math.warp", "Remapping masked image");
496 
497  int const maxCol = destWidth - 1;
498  int const maxRow = destHeight - 1;
499 
500  detail::WarpAtOnePoint<DestImageT, SrcImageT> warpAtOnePoint(srcImage, control, padValue);
501 
502  if (interpLength > 0) {
503  // Use interpolation. Note that 1 produces the same result as no interpolation
504  // but uses this code branch, thus providing an easy way to compare the two branches.
505 
506  // Estimate for number of horizontal interpolation band edges, to reserve memory in vectors
507  int const numColEdges = 2 + ((destWidth - 1) / interpLength);
508 
509  // A list of edge column indices for interpolation bands;
510  // starts at -1, increments by interpLen (except the final interval), and ends at destWidth-1
511  std::vector<int> edgeColList;
512  edgeColList.reserve(numColEdges);
513 
514  // A list of 1/column width for horizontal interpolation bands; the first value is garbage.
515  // The inverse is used for speed because the values are always multiplied.
516  std::vector<double> invWidthList;
517  invWidthList.reserve(numColEdges);
518 
519  // Compute edgeColList and invWidthList
520  edgeColList.push_back(-1);
521  invWidthList.push_back(0.0);
522  for (int prevEndCol = -1; prevEndCol < maxCol; prevEndCol += interpLength) {
523  int endCol = prevEndCol + interpLength;
524  if (endCol > maxCol) {
525  endCol = maxCol;
526  }
527  edgeColList.push_back(endCol);
528  assert(endCol - prevEndCol > 0);
529  invWidthList.push_back(1.0 / static_cast<double>(endCol - prevEndCol));
530  }
531  assert(edgeColList.back() == maxCol);
532 
533  // A list of delta source positions along the edge columns of the horizontal interpolation bands
534  std::vector<lsst::geom::Extent2D> yDeltaSrcPosList(edgeColList.size());
535 
536  // A cache of pixel positions on the source corresponding to the previous or current row
537  // of the destination image.
538  // The first value is for column -1 because the previous source position is used to compute relative
539  // area To simplify the indexing, use an iterator that starts at begin+1, thus: srcPosView =
540  // srcPosList.begin() + 1 srcPosView[col-1] and lower indices are for this row srcPosView[col] and
541  // higher indices are for the previous row
542  std::vector<lsst::geom::Point2D> srcPosList(1 + destWidth);
543  std::vector<lsst::geom::Point2D>::iterator const srcPosView = srcPosList.begin() + 1;
544 
545  std::vector<lsst::geom::Point2D> endColPosList;
546  endColPosList.reserve(numColEdges);
547 
548  // Initialize srcPosList for row -1
549  for (int endCol : edgeColList) {
550  endColPosList.emplace_back(lsst::geom::Point2D(endCol, -1));
551  }
552  auto rightSrcPosList = localDestToParentSrc->applyForward(endColPosList);
553  srcPosView[-1] = rightSrcPosList[0];
554  for (int colBand = 1, endBand = edgeColList.size(); colBand < endBand; ++colBand) {
555  int const prevEndCol = edgeColList[colBand - 1];
556  int const endCol = edgeColList[colBand];
557  lsst::geom::Point2D leftSrcPos = srcPosView[prevEndCol];
558 
559  lsst::geom::Extent2D xDeltaSrcPos =
560  (rightSrcPosList[colBand] - leftSrcPos) * invWidthList[colBand];
561 
562  for (int col = prevEndCol + 1; col <= endCol; ++col) {
563  srcPosView[col] = srcPosView[col - 1] + xDeltaSrcPos;
564  }
565  }
566 
567  int endRow = -1;
568  while (endRow < maxRow) {
569  // Next horizontal interpolation band
570 
571  int prevEndRow = endRow;
572  endRow = prevEndRow + interpLength;
573  if (endRow > maxRow) {
574  endRow = maxRow;
575  }
576  assert(endRow - prevEndRow > 0);
577  double interpInvHeight = 1.0 / static_cast<double>(endRow - prevEndRow);
578 
579  // Set yDeltaSrcPosList for this horizontal interpolation band
580  std::vector<lsst::geom::Point2D> destRowPosList;
581  destRowPosList.reserve(edgeColList.size());
582  for (int endCol : edgeColList) {
583  destRowPosList.emplace_back(lsst::geom::Point2D(endCol, endRow));
584  }
585  auto bottomSrcPosList = localDestToParentSrc->applyForward(destRowPosList);
586  for (int colBand = 0, endBand = edgeColList.size(); colBand < endBand; ++colBand) {
587  int endCol = edgeColList[colBand];
588  yDeltaSrcPosList[colBand] =
589  (bottomSrcPosList[colBand] - srcPosView[endCol]) * interpInvHeight;
590  }
591 
592  for (int row = prevEndRow + 1; row <= endRow; ++row) {
593  typename DestImageT::x_iterator destXIter = destImage.row_begin(row);
594  srcPosView[-1] += yDeltaSrcPosList[0];
595  for (int colBand = 1, endBand = edgeColList.size(); colBand < endBand; ++colBand) {
596  // Next vertical interpolation band
597 
598  int const prevEndCol = edgeColList[colBand - 1];
599  int const endCol = edgeColList[colBand];
600 
601  // Compute xDeltaSrcPos; remember that srcPosView contains
602  // positions for this row in prevEndCol and smaller indices,
603  // and positions for the previous row for larger indices (including endCol)
604  lsst::geom::Point2D leftSrcPos = srcPosView[prevEndCol];
605  lsst::geom::Point2D rightSrcPos = srcPosView[endCol] + yDeltaSrcPosList[colBand];
606  lsst::geom::Extent2D xDeltaSrcPos = (rightSrcPos - leftSrcPos) * invWidthList[colBand];
607 
608  for (int col = prevEndCol + 1; col <= endCol; ++col, ++destXIter) {
609  lsst::geom::Point2D leftSrcPos = srcPosView[col - 1];
610  lsst::geom::Point2D srcPos = leftSrcPos + xDeltaSrcPos;
611  double relativeArea = computeRelativeArea(srcPos, leftSrcPos, srcPosView[col]);
612 
613  srcPosView[col] = srcPos;
614 
615  if (warpAtOnePoint(
616  destXIter, srcPos, relativeArea,
618  ++numGoodPixels;
619  }
620  } // for col
621  } // for col band
622  } // for row
623  } // while next row band
624 
625  } else {
626  // No interpolation
627 
628  // prevSrcPosList = source positions from the previous row; these are used to compute pixel area;
629  // to begin, compute sources positions corresponding to destination row = -1
631  destPosList.reserve(1 + destWidth);
632  for (int col = -1; col < destWidth; ++col) {
633  destPosList.emplace_back(lsst::geom::Point2D(col, -1));
634  }
635  auto prevSrcPosList = localDestToParentSrc->applyForward(destPosList);
636 
637  for (int row = 0; row < destHeight; ++row) {
638  destPosList.clear();
639  for (int col = -1; col < destWidth; ++col) {
640  destPosList.emplace_back(lsst::geom::Point2D(col, row));
641  }
642  auto srcPosList = localDestToParentSrc->applyForward(destPosList);
643 
644  typename DestImageT::x_iterator destXIter = destImage.row_begin(row);
645  for (int col = 0; col < destWidth; ++col, ++destXIter) {
646  // column index = column + 1 because the first entry in srcPosList is for column -1
647  auto srcPos = srcPosList[col + 1];
648  double relativeArea =
649  computeRelativeArea(srcPos, prevSrcPosList[col], prevSrcPosList[col + 1]);
650 
651  if (warpAtOnePoint(destXIter, srcPos, relativeArea,
653  ++numGoodPixels;
654  }
655  } // for col
656  // move points from srcPosList to prevSrcPosList (we don't care about what ends up in srcPosList
657  // because it will be reallocated anyway)
658  swap(srcPosList, prevSrcPosList);
659  } // for row
660  } // if interp
661 
662  return numGoodPixels;
663 }
664 
665 template <typename DestImageT, typename SrcImageT>
666 int warpCenteredImage(DestImageT &destImage, SrcImageT const &srcImage,
667  lsst::geom::LinearTransform const &linearTransform,
668  lsst::geom::Point2D const &centerPosition, WarpingControl const &control,
669  typename DestImageT::SinglePixel padValue) {
670  // force src and dest to be the same size and xy0
671  if ((destImage.getWidth() != srcImage.getWidth()) || (destImage.getHeight() != srcImage.getHeight()) ||
672  (destImage.getXY0() != srcImage.getXY0())) {
673  std::ostringstream errStream;
674  errStream << "src and dest images must have same size and xy0.";
676  }
677 
678  // set the xy0 coords to 0,0 to make life easier
679  SrcImageT srcImageCopy(srcImage, true);
680  srcImageCopy.setXY0(0, 0);
681  destImage.setXY0(0, 0);
682  lsst::geom::Extent2D cLocal =
683  lsst::geom::Extent2D(centerPosition) - lsst::geom::Extent2D(srcImage.getXY0());
684 
685  // for the affine transform, the centerPosition will not only get sheared, but also
686  // moved slightly. So we'll include a translation to move it back by an amount
687  // centerPosition - translatedCenterPosition
688  lsst::geom::AffineTransform affTran(linearTransform, cLocal - linearTransform(cLocal));
690 
691 // now warp
692 #if 0
693  static float t = 0.0;
694  float t_before = 1.0*clock()/CLOCKS_PER_SEC;
695  int n = warpImage(destImage, srcImageCopy, affTran, control, padValue);
696  float t_after = 1.0*clock()/CLOCKS_PER_SEC;
697  float dt = t_after - t_before;
698  t += dt;
699  std::cout <<srcImage.getWidth()<<"x"<<srcImage.getHeight()<<": "<< dt <<" "<< t <<std::endl;
700 #else
701  int n = warpImage(destImage, srcImageCopy, *affineTransform22, control, padValue);
702 #endif
703 
704  // fix the origin and we're done.
705  destImage.setXY0(srcImage.getXY0());
706 
707  return n;
708 }
709 
710 //
711 // Explicit instantiations
712 //
714 // may need to omit default params for EXPOSURE -- original code did that and it worked
715 #define EXPOSURE(PIXTYPE) image::Exposure<PIXTYPE, image::MaskPixel, image::VariancePixel>
716 #define MASKEDIMAGE(PIXTYPE) image::MaskedImage<PIXTYPE, image::MaskPixel, image::VariancePixel>
717 #define IMAGE(PIXTYPE) image::Image<PIXTYPE>
718 #define NL /* */
719 
720 #define INSTANTIATE(DESTIMAGEPIXELT, SRCIMAGEPIXELT) \
721  template int warpCenteredImage( \
722  IMAGE(DESTIMAGEPIXELT) & destImage, IMAGE(SRCIMAGEPIXELT) const &srcImage, \
723  lsst::geom::LinearTransform const &linearTransform, lsst::geom::Point2D const &centerPosition, \
724  WarpingControl const &control, IMAGE(DESTIMAGEPIXELT)::SinglePixel padValue); \
725  NL template int warpCenteredImage( \
726  MASKEDIMAGE(DESTIMAGEPIXELT) & destImage, MASKEDIMAGE(SRCIMAGEPIXELT) const &srcImage, \
727  lsst::geom::LinearTransform const &linearTransform, lsst::geom::Point2D const &centerPosition, \
728  WarpingControl const &control, MASKEDIMAGE(DESTIMAGEPIXELT)::SinglePixel padValue); \
729  NL template int warpImage(IMAGE(DESTIMAGEPIXELT) & destImage, IMAGE(SRCIMAGEPIXELT) const &srcImage, \
730  geom::TransformPoint2ToPoint2 const &srcToDest, WarpingControl const &control, \
731  IMAGE(DESTIMAGEPIXELT)::SinglePixel padValue); \
732  NL template int warpImage(MASKEDIMAGE(DESTIMAGEPIXELT) & destImage, \
733  MASKEDIMAGE(SRCIMAGEPIXELT) const &srcImage, \
734  geom::TransformPoint2ToPoint2 const &srcToDest, WarpingControl const &control, \
735  MASKEDIMAGE(DESTIMAGEPIXELT)::SinglePixel padValue); \
736  NL template int warpImage(IMAGE(DESTIMAGEPIXELT) & destImage, geom::SkyWcs const &destWcs, \
737  IMAGE(SRCIMAGEPIXELT) const &srcImage, geom::SkyWcs const &srcWcs, \
738  WarpingControl const &control, IMAGE(DESTIMAGEPIXELT)::SinglePixel padValue); \
739  NL template int warpImage(MASKEDIMAGE(DESTIMAGEPIXELT) & destImage, geom::SkyWcs const &destWcs, \
740  MASKEDIMAGE(SRCIMAGEPIXELT) const &srcImage, geom::SkyWcs const &srcWcs, \
741  WarpingControl const &control, \
742  MASKEDIMAGE(DESTIMAGEPIXELT)::SinglePixel padValue); \
743  NL template int warpExposure(EXPOSURE(DESTIMAGEPIXELT) & destExposure, \
744  EXPOSURE(SRCIMAGEPIXELT) const &srcExposure, WarpingControl const &control, \
745  EXPOSURE(DESTIMAGEPIXELT)::MaskedImageT::SinglePixel padValue);
746 
747 INSTANTIATE(double, double)
748 INSTANTIATE(double, float)
749 INSTANTIATE(double, int)
750 INSTANTIATE(double, std::uint16_t)
751 INSTANTIATE(float, float)
752 INSTANTIATE(float, int)
754 INSTANTIATE(int, int)
757 } // namespace math
758 
759 template std::shared_ptr<math::LanczosWarpingKernel> table::io::PersistableFacade<
760  math::LanczosWarpingKernel>::dynamicCast(std::shared_ptr<table::io::Persistable> const &);
761 template std::shared_ptr<math::BilinearWarpingKernel> table::io::PersistableFacade<
762  math::BilinearWarpingKernel>::dynamicCast(std::shared_ptr<table::io::Persistable> const &);
763 template std::shared_ptr<math::NearestWarpingKernel> table::io::PersistableFacade<
764  math::NearestWarpingKernel>::dynamicCast(std::shared_ptr<table::io::Persistable> const &);
765 template std::shared_ptr<math::WarpingControl> table::io::PersistableFacade<
766  math::WarpingControl>::dynamicCast(std::shared_ptr<table::io::Persistable> const &);
767 
768 } // namespace afw
769 } // 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:515
Implementation of the Photometric Calibration class.
std::ostream * os
Definition: Schema.cc:557
std::string prefix
Definition: SchemaMapper.cc:72
int y
Definition: SpanSet.cc:48
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
Definition: Persistable.h:48
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.
BilinearWarpingKernel & operator=(const BilinearWarpingKernel &)=delete
void write(OutputArchiveHandle &handle) const override
Write the object to one or more catalogs.
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:212
lsst::geom::Point2I getCtr() const
Return index of kernel's center.
Definition: Kernel.h:234
int getWidth() const
Return the Kernel's width.
Definition: Kernel.h:224
Lanczos warping: accurate but slow and can introduce ringing artifacts.
Definition: warpExposure.h:66
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
void write(OutputArchiveHandle &handle) const override
Write the object to one or more catalogs.
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:180
std::shared_ptr< Kernel > clone() const override
Return a pointer to a deep copy of this kernel.
void write(OutputArchiveHandle &handle) const override
Write the object to one or more catalogs.
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:860
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.
SeparableKernel()
Construct an empty spatially invariant SeparableKernel of size 0x0.
Parameters to control convolution.
Definition: warpExposure.h:276
void setWarpingKernel(SeparableKernel const &warpingKernel)
set the warping kernel
int getInterpLength() const
get the interpolation length (pixels)
Definition: warpExposure.h:329
std::string getPythonModule() const override
Return the fully-qualified Python module that should be imported to guarantee that its factory is reg...
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
std::string getPersistenceName() const override
Return the unique name used to persist this object and look up its factory.
void setMaskWarpingKernel(SeparableKernel const &maskWarpingKernel)
set the mask warping kernel
std::shared_ptr< SeparableKernel > getWarpingKernel() const
get the warping kernel
void write(OutputArchiveHandle &handle) const override
Write the object to one or more catalogs.
std::shared_ptr< SeparableKernel > getMaskWarpingKernel() const
get the mask warping kernel
bool isPersistable() const noexcept override
Return true if this particular object can be persisted using afw::table::io.
A functor that computes one warped pixel.
std::shared_ptr< RecordT > addNew()
Create a new record, add it to the end of the catalog, and return a pointer to it.
Definition: Catalog.h:490
Defines the fields and offsets for a table.
Definition: Schema.h:51
An object passed to Persistable::write to allow it to persist itself.
void saveCatalog(BaseCatalog const &catalog)
Save a catalog in the archive.
void saveEmpty()
Indicate that the object being persisted has no state, and hence will never call makeCatalog() or sav...
BaseCatalog makeCatalog(Schema const &schema)
Return a new, empty catalog with the given schema.
int put(Persistable const *obj, bool permissive=false)
Save an object to the archive and return a unique ID that can be used to retrieve it from an InputArc...
PersistableFactory(std::string const &name)
Constructor for the factory.
Definition: Persistable.cc:74
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:146
Transform< Point2Endpoint, Point2Endpoint > TransformPoint2ToPoint2
Definition: Transform.h:300
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
std::int32_t MaskPixel
default type for Masks and MaskedImage Masks
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:692
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 regex_match(T... args)
T reserve(T... args)
T size(T... args)
afw::table::Key< std::string > warpingKernelName
Definition: CoaddPsf.cc:350
int row
Definition: CR.cc:145
int col
Definition: CR.cc:144
T stoi(T... args)
T str(T... args)
typename ImageT::image_category image_category
Definition: ImageBase.h:67
T swap(T... args)
table::Key< int > warpingKernelIndex
table::Schema schema
std::shared_ptr< table::io::Persistable > read(table::io::InputArchive const &archive, table::io::CatalogVector const &catalogs) const override
Definition: warpExposure.cc:0
table::Key< int > interpLength
table::Key< image::MaskPixel > growFullMask
table::Key< int > order
table::Key< int > cacheSize
table::Key< int > maskKernelIndex
table::Key< table::Flag > hasMaskKernel