LSST Applications g0b6bd0c080+a72a5dd7e6,g1182afd7b4+2a019aa3bb,g17e5ecfddb+2b8207f7de,g1d67935e3f+06cf436103,g38293774b4+ac198e9f13,g396055baef+6a2097e274,g3b44f30a73+6611e0205b,g480783c3b1+98f8679e14,g48ccf36440+89c08d0516,g4b93dc025c+98f8679e14,g5c4744a4d9+a302e8c7f0,g613e996a0d+e1c447f2e0,g6c8d09e9e7+25247a063c,g7271f0639c+98f8679e14,g7a9cd813b8+124095ede6,g9d27549199+a302e8c7f0,ga1cf026fa3+ac198e9f13,ga32aa97882+7403ac30ac,ga786bb30fb+7a139211af,gaa63f70f4e+9994eb9896,gabf319e997+ade567573c,gba47b54d5d+94dc90c3ea,gbec6a3398f+06cf436103,gc6308e37c7+07dd123edb,gc655b1545f+ade567573c,gcc9029db3c+ab229f5caf,gd01420fc67+06cf436103,gd877ba84e5+06cf436103,gdb4cecd868+6f279b5b48,ge2d134c3d5+cc4dbb2e3f,ge448b5faa6+86d1ceac1d,gecc7e12556+98f8679e14,gf3ee170dca+25247a063c,gf4ac96e456+ade567573c,gf9f5ea5b4d+ac198e9f13,gff490e6085+8c2580be5c,w.2022.27
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"
51#include "lsst/afw/table/io/Persistable.cc" // Needed for PersistableFacade::dynamicCast
52
54
55using std::swap;
56
57namespace lsst {
58namespace afw {
59namespace 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//
73static 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
101int LanczosWarpingKernel::getOrder() const { return this->getWidth() / 2; }
102
103void 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
128void 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
149void 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
161namespace {
162
163struct 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
177private:
178 LanczosKernelPersistenceHelper()
179 : schema(), order(schema.addField<int>("order", "order of Lanczos function")) {}
180};
181
182class : 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
196template <class T>
197class 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
207DefaultPersistableFactory<BilinearWarpingKernel> bilinearFactory("BilinearWarpingKernel");
208DefaultPersistableFactory<NearestWarpingKernel> nearestFactory("NearestWarpingKernel");
209
210} // namespace
211
213 auto const &keys = LanczosKernelPersistenceHelper::get();
214 table::BaseCatalog catalog = handle.makeCatalog(keys.schema);
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
271void 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
285void 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
299namespace {
300
301struct 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
320private:
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")),
333 "growFullMask", "bits to grow to full width of image/variance kernel")) {}
334};
335
336std::string _getWarpingControlPersistenceName() { return "WarpingControl"; }
337
338class : 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
366std::string WarpingControl::getPersistenceName() const { return _getWarpingControlPersistenceName(); }
367
368std::string WarpingControl::getPythonModule() const { return "lsst.afw.math"; }
369
370bool 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);
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
392template <typename DestExposureT, typename SrcExposureT>
393int 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.setFilter(srcExposure.getFilter());
407 destExposure.getInfo()->setVisitInfo(srcExposure.getInfo()->getVisitInfo());
408 return warpImage(mi, *destExposure.getWcs(), srcExposure.getMaskedImage(), *srcExposure.getWcs(), control,
409 padValue);
410}
411
412namespace {
413
414inline 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
426inline double computeRelativeArea(
427 lsst::geom::Point2D const &srcPos,
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
440template <typename DestImageT, typename SrcImageT>
441int 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
448template <typename DestImageT, typename SrcImageT>
449int 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
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
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) {
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
665template <typename DestImageT, typename SrcImageT>
666int 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
747INSTANTIATE(double, double)
748INSTANTIATE(double, float)
749INSTANTIATE(double, int)
751INSTANTIATE(float, float)
752INSTANTIATE(float, int)
754INSTANTIATE(int, int)
757} // namespace math
758
759template std::shared_ptr<math::LanczosWarpingKernel> table::io::PersistableFacade<
760 math::LanczosWarpingKernel>::dynamicCast(std::shared_ptr<table::io::Persistable> const &);
761template std::shared_ptr<math::BilinearWarpingKernel> table::io::PersistableFacade<
762 math::BilinearWarpingKernel>::dynamicCast(std::shared_ptr<table::io::Persistable> const &);
763template std::shared_ptr<math::NearestWarpingKernel> table::io::PersistableFacade<
764 math::NearestWarpingKernel>::dynamicCast(std::shared_ptr<table::io::Persistable> const &);
765template 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:733
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:354
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