LSST Applications g04c3c9f7ca+2075667efa,g1e125bf412+5f448d5fcf,g2079a07aa2+3e9fd84d81,g2305ad1205+b635cf1488,g2bbee38e9b+6c6beb4891,g337abbeb29+6c6beb4891,g33d1c0ed96+6c6beb4891,g3a166c0a6a+6c6beb4891,g50ff169b8f+96c6868917,g52b1c1532d+585e252eca,g591dd9f2cf+42f171e1e6,g5c3423f6d4+d536b04327,g607f77f49a+d536b04327,g6f43f06aed+ca1339dc19,g858d7b2824+d536b04327,g8ee334c5b4+d7f9608c2f,g9963eaa53e+b3dc1655d3,g998f4353bf+d536b04327,g99cad8db69+8ef2408349,g9ddcbc5298+9a081db1e4,ga1e77700b3+2cbb763275,gadfd92a7e4+aec2f3b930,gae0086650b+585e252eca,gb0e22166c9+0e73c8378f,gb3b7280ab2+cb5fdb229e,gbb8dafda3b+a327199e22,gc120e1dc64+88074880ea,gc28159a63d+6c6beb4891,gcdd4ae20e8+bd241b2308,gcde1bda545+903e937d91,gcf0d15dbbd+bd241b2308,gdaeeff99f8+f9a426f77a,gddc38dedce+585e252eca,ge79ae78c31+6c6beb4891,gfbcc870c63+b310236976,w.2024.23
LSST Data Management Base Package
Loading...
Searching...
No Matches
Namespaces | Functions
warpExposure.cc File Reference
#include <cassert>
#include <cmath>
#include <limits>
#include <memory>
#include <sstream>
#include <string>
#include <vector>
#include <utility>
#include <regex>
#include "astshim.h"
#include "lsst/log/Log.h"
#include "lsst/pex/exceptions.h"
#include "lsst/geom.h"
#include "lsst/afw/math/warpExposure.h"
#include "lsst/afw/geom.h"
#include "lsst/afw/math/Kernel.h"
#include "lsst/afw/image/PhotoCalib.h"
#include "lsst/afw/math/detail/WarpAtOnePoint.h"
#include "lsst/afw/table/io/InputArchive.h"
#include "lsst/afw/table/io/OutputArchive.h"
#include "lsst/afw/table/io/CatalogVector.h"
#include "lsst/afw/table/io/Persistable.cc"

Go to the source code of this file.

Namespaces

namespace  lsst
 
namespace  lsst::afw
 
namespace  lsst::afw::math
 

Functions

std::shared_ptr< table::io::Persistable > read (table::io::InputArchive const &archive, table::io::CatalogVector const &catalogs) const override
 
std::shared_ptr< SeparableKernellsst::afw::math::makeWarpingKernel (std::string name)
 Return a warping kernel given its name.
 
template<typename DestExposureT , typename SrcExposureT >
int lsst::afw::math::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.
 
template<typename DestImageT , typename SrcImageT >
int lsst::afw::math::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.
 
template<typename DestImageT , typename SrcImageT >
int lsst::afw::math::warpImage (DestImageT &destImage, SrcImageT const &srcImage, geom::TransformPoint2ToPoint2 const &srcToDest, WarpingControl const &control, typename DestImageT::SinglePixel padValue=lsst::afw::math::edgePixel< DestImageT >(typename lsst::afw::image::detail::image_traits< DestImageT >::image_category()))
 A variant of warpImage that uses a Transform<Point2Endpoint, Point2Endpoint> instead of a pair of WCS to describe the transformation.
 
template<typename DestImageT , typename SrcImageT >
int lsst::afw::math::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.
 

Function Documentation

◆ read()

std::shared_ptr< table::io::Persistable > read ( table::io::InputArchive const & archive,
table::io::CatalogVector const & catalogs ) const
inlineoverride

Definition at line 0 of file warpExposure.cc.

9 : you can redistribute it and/or modify
10 * it under the terms of the GNU General Public License as published by
11 * the Free Software Foundation, either version 3 of the License, or
12 * (at your option) any later version.
13 *
14 * This program is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 * GNU General Public License for more details.
18 *
19 * You should have received a copy of the LSST License Statement and
20 * the GNU General Public License along with this program. If not,
21 * see <http://www.lsstcorp.org/LegalNotices/>.
22 */
23/*
24 * Support for warping an %image to a new Wcs.
25 */
26
27#include <cassert>
28#include <cmath>
29#include <limits>
30#include <memory>
31#include <sstream>
32#include <string>
33#include <vector>
34#include <utility>
35#include <regex>
36
37#include "astshim.h"
38
39#include "lsst/log/Log.h"
40#include "lsst/pex/exceptions.h"
41#include "lsst/geom.h"
43#include "lsst/afw/geom.h"
50#include "lsst/afw/table/io/Persistable.cc" // Needed for PersistableFacade::dynamicCast
51
53
54using std::swap;
55
56namespace lsst {
57namespace afw {
58namespace math {
59
60//
61// A helper function for the warping kernels which provides error-checking:
62// the warping kernels are designed to work in two cases
63// 0 < x < 1 and ctrX=(size-1)/2
64// -1 < x < 0 and ctrX=(size+1)/2
65// (and analogously for y). Note that to get the second case, Kernel::setCtr(1, y) must be
66// called before calling Kernel::setKernelParameter(). [see afw::math::offsetImage() for
67// an example]
68//
69// FIXME eventually the 3 warping kernels will inherit from a common base class WarpingKernel
70// and this routine can be eliminated by putting the code in WarpingKernel::setKernelParameter()
71//
72static inline void checkWarpingKernelParameter(const SeparableKernel *p, unsigned int ind, double value) {
73 if (ind > 1) {
75 "bad ind argument in WarpingKernel::setKernelParameter()");
76 }
77 int ctr = p->getCtr()[ind];
78 int size = p->getDimensions()[ind];
79
80 if (ctr == (size - 1) / 2) {
81 if (value < -1e-6 || value > 1 + 1e-6) {
83 "bad coordinate in WarpingKernel::setKernelParameter()");
84 }
85 } else if (ctr == (size + 1) / 2) {
86 if (value < -1 - 1e-6 || value > 1e-6) {
88 "bad coordinate in WarpingKernel::setKernelParameter()");
89 }
90 } else {
92 "bad ctr value in WarpingKernel::setKernelParameter()");
93 }
94}
95
98}
99
100int LanczosWarpingKernel::getOrder() const { return this->getWidth() / 2; }
101
102void LanczosWarpingKernel::setKernelParameter(unsigned int ind, double value) const {
103 checkWarpingKernelParameter(this, ind, value);
105}
106
109}
110
112 //
113 // this->_params[0] = value of x where we want to interpolate the function
114 // x = integer value of x where we evaluate the function in the interpolation
115 //
116 // The following weird-looking expression has no if/else statements, is roundoff-tolerant,
117 // and works in the following two cases:
118 // 0 < this->_params[0] < 1, x \in {0,1}
119 // -1 < this->_params[0] < 0, x \in {-1,0}
120 //
121 // The checks in BilinearWarpingKernel::setKernelParameter() ensure that one of these
122 // conditions is satisfied
123 //
124 return 0.5 + (1.0 - (2.0 * fabs(this->_params[0]))) * (0.5 - fabs(x));
125}
126
127void BilinearWarpingKernel::setKernelParameter(unsigned int ind, double value) const {
128 checkWarpingKernelParameter(this, ind, value);
130}
131
134 os << "_BilinearFunction1: ";
135 os << Function1<Kernel::Pixel>::toString(prefix);
136 return os.str();
137}
138
140 return std::make_shared<NearestWarpingKernel>();
141}
142
144 // this expression is faster than using conditionals, but offers no sanity checking
145 return static_cast<double>((fabs(this->_params[0]) < 0.5) == (fabs(x) < 0.5));
146}
147
148void NearestWarpingKernel::setKernelParameter(unsigned int ind, double value) const {
149 checkWarpingKernelParameter(this, ind, value);
151}
152
155 os << "_NearestFunction1: ";
156 os << Function1<Kernel::Pixel>::toString(prefix);
157 return os.str();
158}
159
160namespace {
161
162struct LanczosKernelPersistenceHelper {
163 table::Schema schema;
164 table::Key<int> order;
165
166 static LanczosKernelPersistenceHelper const &get() {
167 static LanczosKernelPersistenceHelper const instance;
168 return instance;
169 }
170
171 LanczosKernelPersistenceHelper(LanczosKernelPersistenceHelper const &) = delete;
172 LanczosKernelPersistenceHelper(LanczosKernelPersistenceHelper &&) = delete;
173 LanczosKernelPersistenceHelper &operator=(LanczosKernelPersistenceHelper const &) = delete;
174 LanczosKernelPersistenceHelper &operator=(LanczosKernelPersistenceHelper &&) = delete;
175
176private:
177 LanczosKernelPersistenceHelper()
178 : schema(), order(schema.addField<int>("order", "order of Lanczos function")) {}
179};
180
181class : public table::io::PersistableFactory {
182 std::shared_ptr<table::io::Persistable> read(table::io::InputArchive const &archive,
183 table::io::CatalogVector const &catalogs) const override {
184 auto const &keys = LanczosKernelPersistenceHelper::get();
185 LSST_ARCHIVE_ASSERT(catalogs.size() == 1u);
186 LSST_ARCHIVE_ASSERT(catalogs.front().size() == 1u);
187 afw::table::BaseRecord const &record = catalogs.front().front();
188 LSST_ARCHIVE_ASSERT(record.getSchema() == keys.schema);
189 return std::make_shared<LanczosWarpingKernel>(record.get(keys.order));
190 }
191
193} lanczosFactory("LanczosWarpingKernel");
194
195template <class T>
196class DefaultPersistableFactory : public table::io::PersistableFactory {
197 std::shared_ptr<table::io::Persistable> read(table::io::InputArchive const &archive,
198 table::io::CatalogVector const &catalogs) const override {
199 LSST_ARCHIVE_ASSERT(catalogs.empty());
200 return std::make_shared<T>();
201 }
202
203 using table::io::PersistableFactory::PersistableFactory;
204};
205
206DefaultPersistableFactory<BilinearWarpingKernel> bilinearFactory("BilinearWarpingKernel");
207DefaultPersistableFactory<NearestWarpingKernel> nearestFactory("NearestWarpingKernel");
208
209} // namespace
210
212 auto const &keys = LanczosKernelPersistenceHelper::get();
213 table::BaseCatalog catalog = handle.makeCatalog(keys.schema);
214 std::shared_ptr<table::BaseRecord> record = catalog.addNew();
215 record->set(keys.order, getOrder());
216 handle.saveCatalog(catalog);
217}
218
219void BilinearWarpingKernel::write(OutputArchiveHandle &handle) const { handle.saveEmpty(); }
220
221void NearestWarpingKernel::write(OutputArchiveHandle &handle) const { handle.saveEmpty(); }
222
224 using KernelPtr = std::shared_ptr<SeparableKernel>;
225 std::smatch matches;
226 static const std::regex LanczosRE("lanczos(\\d+)");
227 if (name == "bilinear") {
228 return KernelPtr(new BilinearWarpingKernel());
229 } else if (std::regex_match(name, matches, LanczosRE)) {
230 std::string orderStr(matches[1].first, matches[1].second);
231 int order = std::stoi(orderStr);
232 return KernelPtr(new LanczosWarpingKernel(order));
233 } else if (name == "nearest") {
234 return KernelPtr(new NearestWarpingKernel());
235 } else {
236 throw LSST_EXCEPT(pexExcept::InvalidParameterError, "unknown warping kernel name: \"" + name + "\"");
237 }
238}
239
241 if (_warpingKernelPtr->getCacheSize() != _cacheSize) {
242 _warpingKernelPtr->computeCache(_cacheSize);
243 }
244 return _warpingKernelPtr;
245};
246
249 setWarpingKernel(*warpingKernelPtr);
250}
251
252void WarpingControl::setWarpingKernel(SeparableKernel const &warpingKernel) {
253 if (_maskWarpingKernelPtr) {
254 _testWarpingKernels(warpingKernel, *_maskWarpingKernelPtr);
255 }
256 std::shared_ptr<SeparableKernel> warpingKernelPtr(
257 std::static_pointer_cast<SeparableKernel>(warpingKernel.clone()));
258 _warpingKernelPtr = warpingKernelPtr;
259}
260
262 if (_maskWarpingKernelPtr) { // lazily update kernel cache
263 if (_maskWarpingKernelPtr->getCacheSize() != _cacheSize) {
264 _maskWarpingKernelPtr->computeCache(_cacheSize);
265 }
266 }
267 return _maskWarpingKernelPtr;
268}
269
270void WarpingControl::setMaskWarpingKernelName(std::string const &maskWarpingKernelName) {
271 if (!maskWarpingKernelName.empty()) {
272 std::shared_ptr<SeparableKernel> maskWarpingKernelPtr(makeWarpingKernel(maskWarpingKernelName));
273 setMaskWarpingKernel(*maskWarpingKernelPtr);
274 } else {
275 _maskWarpingKernelPtr.reset();
276 }
277}
278
279void WarpingControl::setMaskWarpingKernel(SeparableKernel const &maskWarpingKernel) {
280 _testWarpingKernels(*_warpingKernelPtr, maskWarpingKernel);
281 _maskWarpingKernelPtr = std::static_pointer_cast<SeparableKernel>(maskWarpingKernel.clone());
282}
283
284void WarpingControl::_testWarpingKernels(SeparableKernel const &warpingKernel,
285 SeparableKernel const &maskWarpingKernel) const {
286 lsst::geom::Box2I kernelBBox =
287 lsst::geom::Box2I(lsst::geom::Point2I(0, 0) - lsst::geom::Extent2I(warpingKernel.getCtr()),
288 warpingKernel.getDimensions());
289 lsst::geom::Box2I maskKernelBBox =
290 lsst::geom::Box2I(lsst::geom::Point2I(0, 0) - lsst::geom::Extent2I(maskWarpingKernel.getCtr()),
291 maskWarpingKernel.getDimensions());
292 if (!kernelBBox.contains(maskKernelBBox)) {
293 throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
294 "warping kernel is smaller than mask warping kernel");
295 }
296}
297
298namespace {
299
300struct WarpingControlPersistenceHelper {
301 table::Schema schema;
302 table::Key<int> warpingKernelIndex;
303 table::Key<table::Flag> hasMaskKernel;
304 table::Key<int> maskKernelIndex; // undefined if !hasMaskKernel
305 table::Key<int> cacheSize;
306 table::Key<int> interpLength;
307 table::Key<image::MaskPixel> growFullMask;
308
309 static WarpingControlPersistenceHelper const &get() {
310 static WarpingControlPersistenceHelper const instance;
311 return instance;
312 }
313
314 WarpingControlPersistenceHelper(WarpingControlPersistenceHelper const &) = delete;
315 WarpingControlPersistenceHelper(WarpingControlPersistenceHelper &&) = delete;
316 WarpingControlPersistenceHelper &operator=(WarpingControlPersistenceHelper const &) = delete;
317 WarpingControlPersistenceHelper &operator=(WarpingControlPersistenceHelper &&) = delete;
318
319private:
320 WarpingControlPersistenceHelper()
321 : schema(),
323 schema.addField<int>("warpingKernelIndex", "archive ID of nested warping kernel")),
324 hasMaskKernel(schema.addField<table::Flag>("hasMaskKernel", "whether a mask kernel is stored")),
325 maskKernelIndex(schema.addField<int>("maskKernelIndex",
326 "archive ID of nested mask kernel. "
327 "Valid only if hasMaskKernel")),
328 cacheSize(schema.addField<int>("cacheSize", "Cache size for warping kernel(s)")),
329 interpLength(schema.addField<int>("interpLength",
330 "Distance over which WCS can be linearly interpolated")),
331 growFullMask(schema.addField<image::MaskPixel>(
332 "growFullMask", "bits to grow to full width of image/variance kernel")) {}
333};
334
335std::string _getWarpingControlPersistenceName() { return "WarpingControl"; }
336
337class : public table::io::PersistableFactory {
338 std::shared_ptr<table::io::Persistable> read(table::io::InputArchive const &archive,
339 table::io::CatalogVector const &catalogs) const override {
340 auto const &keys = WarpingControlPersistenceHelper::get();
341 LSST_ARCHIVE_ASSERT(catalogs.size() == 1u);
342 LSST_ARCHIVE_ASSERT(catalogs.front().size() == 1u);
343 afw::table::BaseRecord const &record = catalogs.front().front();
344 LSST_ARCHIVE_ASSERT(record.getSchema() == keys.schema);
345
346 // Warping kernels are dummy values, set true kernels later
347 auto control = std::make_shared<WarpingControl>("bilinear", "", record.get(keys.cacheSize),
348 record.get(keys.interpLength),
349 record.get(keys.growFullMask));
350 // archive.get returns a shared_ptr, so this code depends on the
351 // undocumented fact that setWarpingKernel and setMaskWarpingKernel
352 // make defensive copies.
353 control->setWarpingKernel(*archive.get<SeparableKernel>(record.get(keys.warpingKernelIndex)));
354 if (record.get(keys.hasMaskKernel)) {
355 control->setMaskWarpingKernel(*archive.get<SeparableKernel>(record.get(keys.maskKernelIndex)));
356 }
357 return control;
358 }
359
361} controlFactory(_getWarpingControlPersistenceName());
362
363} // namespace
364
365std::string WarpingControl::getPersistenceName() const { return _getWarpingControlPersistenceName(); }
366
367std::string WarpingControl::getPythonModule() const { return "lsst.afw.math"; }
368
369bool WarpingControl::isPersistable() const noexcept {
370 return _warpingKernelPtr->isPersistable() &&
371 (!hasMaskWarpingKernel() || _maskWarpingKernelPtr->isPersistable());
372}
373
374void WarpingControl::write(OutputArchiveHandle &handle) const {
375 auto const &keys = WarpingControlPersistenceHelper::get();
376 table::BaseCatalog catalog = handle.makeCatalog(keys.schema);
377 std::shared_ptr<table::BaseRecord> record = catalog.addNew();
378
379 record->set(keys.warpingKernelIndex, handle.put(_warpingKernelPtr));
380 record->set(keys.hasMaskKernel, hasMaskWarpingKernel());
381 if (hasMaskWarpingKernel()) {
382 record->set(keys.maskKernelIndex, handle.put(_maskWarpingKernelPtr));
383 }
384 record->set(keys.cacheSize, _cacheSize);
385 record->set(keys.interpLength, _interpLength);
386 record->set(keys.growFullMask, _growFullMask);
387
388 handle.saveCatalog(catalog);
389}
390
391template <typename DestExposureT, typename SrcExposureT>
392int warpExposure(DestExposureT &destExposure, SrcExposureT const &srcExposure, WarpingControl const &control,
393 typename DestExposureT::MaskedImageT::SinglePixel padValue) {
394 if (!destExposure.hasWcs()) {
395 throw LSST_EXCEPT(pexExcept::InvalidParameterError, "destExposure has no Wcs");
396 }
397 if (!srcExposure.hasWcs()) {
398 throw LSST_EXCEPT(pexExcept::InvalidParameterError, "srcExposure has no Wcs");
399 }
400 typename DestExposureT::MaskedImageT mi = destExposure.getMaskedImage();
401 if (srcExposure.getInfo()->hasId()) {
402 destExposure.getInfo()->setId(srcExposure.getInfo()->getId());
403 }
404 destExposure.setPhotoCalib(srcExposure.getPhotoCalib());
405 destExposure.setFilter(srcExposure.getFilter());
406 destExposure.getInfo()->setVisitInfo(srcExposure.getInfo()->getVisitInfo());
407 return warpImage(mi, *destExposure.getWcs(), srcExposure.getMaskedImage(), *srcExposure.getWcs(), control,
408 padValue);
409}
410
411namespace {
412
413inline lsst::geom::Point2D computeSrcPos(
414 int destCol,
415 int destRow,
416 lsst::geom::Point2D const &destXY0,
417 geom::SkyWcs const &destWcs,
418 geom::SkyWcs const &srcWcs)
419{
420 lsst::geom::Point2D const destPix(image::indexToPosition(destCol + destXY0[0]),
421 image::indexToPosition(destRow + destXY0[1]));
422 return srcWcs.skyToPixel(destWcs.pixelToSky(destPix));
423}
424
425inline double computeRelativeArea(
426 lsst::geom::Point2D const &srcPos,
428 &leftSrcPos,
429 lsst::geom::Point2D const &upSrcPos)
430{
431 lsst::geom::Extent2D dSrcA = srcPos - leftSrcPos;
432 lsst::geom::Extent2D dSrcB = srcPos - upSrcPos;
433
434 return std::abs(dSrcA.getX() * dSrcB.getY() - dSrcA.getY() * dSrcB.getX());
435}
436
437} // namespace
438
439template <typename DestImageT, typename SrcImageT>
440int warpImage(DestImageT &destImage, geom::SkyWcs const &destWcs, SrcImageT const &srcImage,
441 geom::SkyWcs const &srcWcs, WarpingControl const &control,
442 typename DestImageT::SinglePixel padValue) {
443 auto srcToDest = geom::makeWcsPairTransform(srcWcs, destWcs);
444 return warpImage(destImage, srcImage, *srcToDest, control, padValue);
445}
446
447template <typename DestImageT, typename SrcImageT>
448int warpImage(DestImageT &destImage, SrcImageT const &srcImage,
449 geom::TransformPoint2ToPoint2 const &srcToDest, WarpingControl const &control,
450 typename DestImageT::SinglePixel padValue) {
451 if (imagesOverlap(destImage, srcImage)) {
452 throw LSST_EXCEPT(pexExcept::InvalidParameterError, "destImage overlaps srcImage; cannot warp");
453 }
454 if (destImage.getBBox(image::LOCAL).isEmpty()) {
455 return 0;
456 }
457 // if src image is too small then don't try to warp
458 std::shared_ptr<SeparableKernel> warpingKernelPtr = control.getWarpingKernel();
459 try {
460 warpingKernelPtr->shrinkBBox(srcImage.getBBox(image::LOCAL));
462 for (int y = 0, height = destImage.getHeight(); y < height; ++y) {
463 for (typename DestImageT::x_iterator destPtr = destImage.row_begin(y), end = destImage.row_end(y);
464 destPtr != end; ++destPtr) {
465 *destPtr = padValue;
466 }
467 }
468 return 0;
469 }
470 int interpLength = control.getInterpLength();
471
472 std::shared_ptr<LanczosWarpingKernel const> const lanczosKernelPtr =
473 std::dynamic_pointer_cast<LanczosWarpingKernel>(warpingKernelPtr);
474
475 int numGoodPixels = 0;
476
477 // compute a transform from local destination pixels to parent source pixels
478 auto const parentDestToParentSrc = srcToDest.inverted();
479 std::vector<double> const localDestToParentDestVec = {static_cast<double>(destImage.getX0()),
480 static_cast<double>(destImage.getY0())};
481 auto const localDestToParentDest = geom::TransformPoint2ToPoint2(ast::ShiftMap(localDestToParentDestVec));
482 auto const localDestToParentSrc = localDestToParentDest.then(*parentDestToParentSrc);
483
484 // Get the source MaskedImage and a pixel accessor to it.
485 int const srcWidth = srcImage.getWidth();
486 int const srcHeight = srcImage.getHeight();
487 LOGL_DEBUG("TRACE2.lsst.afw.math.warp", "source image width=%d; height=%d", srcWidth, srcHeight);
488
489 int const destWidth = destImage.getWidth();
490 int const destHeight = destImage.getHeight();
491 LOGL_DEBUG("TRACE2.lsst.afw.math.warp", "remap image width=%d; height=%d", destWidth, destHeight);
492
493 // Set each pixel of destExposure's MaskedImage
494 LOGL_DEBUG("TRACE3.lsst.afw.math.warp", "Remapping masked image");
495
496 int const maxCol = destWidth - 1;
497 int const maxRow = destHeight - 1;
498
499 detail::WarpAtOnePoint<DestImageT, SrcImageT> warpAtOnePoint(srcImage, control, padValue);
500
501 if (interpLength > 0) {
502 // Use interpolation. Note that 1 produces the same result as no interpolation
503 // but uses this code branch, thus providing an easy way to compare the two branches.
504
505 // Estimate for number of horizontal interpolation band edges, to reserve memory in vectors
506 int const numColEdges = 2 + ((destWidth - 1) / interpLength);
507
508 // A list of edge column indices for interpolation bands;
509 // starts at -1, increments by interpLen (except the final interval), and ends at destWidth-1
510 std::vector<int> edgeColList;
511 edgeColList.reserve(numColEdges);
512
513 // A list of 1/column width for horizontal interpolation bands; the first value is garbage.
514 // The inverse is used for speed because the values are always multiplied.
515 std::vector<double> invWidthList;
516 invWidthList.reserve(numColEdges);
517
518 // Compute edgeColList and invWidthList
519 edgeColList.push_back(-1);
520 invWidthList.push_back(0.0);
521 for (int prevEndCol = -1; prevEndCol < maxCol; prevEndCol += interpLength) {
522 int endCol = prevEndCol + interpLength;
523 if (endCol > maxCol) {
524 endCol = maxCol;
525 }
526 edgeColList.push_back(endCol);
527 assert(endCol - prevEndCol > 0);
528 invWidthList.push_back(1.0 / static_cast<double>(endCol - prevEndCol));
529 }
530 assert(edgeColList.back() == maxCol);
531
532 // A list of delta source positions along the edge columns of the horizontal interpolation bands
533 std::vector<lsst::geom::Extent2D> yDeltaSrcPosList(edgeColList.size());
534
535 // A cache of pixel positions on the source corresponding to the previous or current row
536 // of the destination image.
537 // The first value is for column -1 because the previous source position is used to compute relative
538 // area To simplify the indexing, use an iterator that starts at begin+1, thus: srcPosView =
539 // srcPosList.begin() + 1 srcPosView[col-1] and lower indices are for this row srcPosView[col] and
540 // higher indices are for the previous row
541 std::vector<lsst::geom::Point2D> srcPosList(1 + destWidth);
542 std::vector<lsst::geom::Point2D>::iterator const srcPosView = srcPosList.begin() + 1;
543
545 endColPosList.reserve(numColEdges);
546
547 // Initialize srcPosList for row -1
548 for (int endCol : edgeColList) {
549 endColPosList.emplace_back(lsst::geom::Point2D(endCol, -1));
550 }
551 auto rightSrcPosList = localDestToParentSrc->applyForward(endColPosList);
552 srcPosView[-1] = rightSrcPosList[0];
553 for (int colBand = 1, endBand = edgeColList.size(); colBand < endBand; ++colBand) {
554 int const prevEndCol = edgeColList[colBand - 1];
555 int const endCol = edgeColList[colBand];
556 lsst::geom::Point2D leftSrcPos = srcPosView[prevEndCol];
557
558 lsst::geom::Extent2D xDeltaSrcPos =
559 (rightSrcPosList[colBand] - leftSrcPos) * invWidthList[colBand];
560
561 for (int col = prevEndCol + 1; col <= endCol; ++col) {
562 srcPosView[col] = srcPosView[col - 1] + xDeltaSrcPos;
563 }
564 }
565
566 int endRow = -1;
567 while (endRow < maxRow) {
568 // Next horizontal interpolation band
569
570 int prevEndRow = endRow;
571 endRow = prevEndRow + interpLength;
572 if (endRow > maxRow) {
573 endRow = maxRow;
574 }
575 assert(endRow - prevEndRow > 0);
576 double interpInvHeight = 1.0 / static_cast<double>(endRow - prevEndRow);
577
578 // Set yDeltaSrcPosList for this horizontal interpolation band
580 destRowPosList.reserve(edgeColList.size());
581 for (int endCol : edgeColList) {
582 destRowPosList.emplace_back(lsst::geom::Point2D(endCol, endRow));
583 }
584 auto bottomSrcPosList = localDestToParentSrc->applyForward(destRowPosList);
585 for (int colBand = 0, endBand = edgeColList.size(); colBand < endBand; ++colBand) {
586 int endCol = edgeColList[colBand];
587 yDeltaSrcPosList[colBand] =
588 (bottomSrcPosList[colBand] - srcPosView[endCol]) * interpInvHeight;
589 }
590
591 for (int row = prevEndRow + 1; row <= endRow; ++row) {
592 typename DestImageT::x_iterator destXIter = destImage.row_begin(row);
593 srcPosView[-1] += yDeltaSrcPosList[0];
594 for (int colBand = 1, endBand = edgeColList.size(); colBand < endBand; ++colBand) {
595 // Next vertical interpolation band
596
597 int const prevEndCol = edgeColList[colBand - 1];
598 int const endCol = edgeColList[colBand];
599
600 // Compute xDeltaSrcPos; remember that srcPosView contains
601 // positions for this row in prevEndCol and smaller indices,
602 // and positions for the previous row for larger indices (including endCol)
603 lsst::geom::Point2D leftSrcPos = srcPosView[prevEndCol];
604 lsst::geom::Point2D rightSrcPos = srcPosView[endCol] + yDeltaSrcPosList[colBand];
605 lsst::geom::Extent2D xDeltaSrcPos = (rightSrcPos - leftSrcPos) * invWidthList[colBand];
606
607 for (int col = prevEndCol + 1; col <= endCol; ++col, ++destXIter) {
608 lsst::geom::Point2D leftSrcPos = srcPosView[col - 1];
609 lsst::geom::Point2D srcPos = leftSrcPos + xDeltaSrcPos;
610 double relativeArea = computeRelativeArea(srcPos, leftSrcPos, srcPosView[col]);
611
612 srcPosView[col] = srcPos;
613
614 if (warpAtOnePoint(
615 destXIter, srcPos, relativeArea,
617 ++numGoodPixels;
618 }
619 } // for col
620 } // for col band
621 } // for row
622 } // while next row band
623
624 } else {
625 // No interpolation
626
627 // prevSrcPosList = source positions from the previous row; these are used to compute pixel area;
628 // to begin, compute sources positions corresponding to destination row = -1
630 destPosList.reserve(1 + destWidth);
631 for (int col = -1; col < destWidth; ++col) {
632 destPosList.emplace_back(lsst::geom::Point2D(col, -1));
633 }
634 auto prevSrcPosList = localDestToParentSrc->applyForward(destPosList);
635
636 for (int row = 0; row < destHeight; ++row) {
637 destPosList.clear();
638 for (int col = -1; col < destWidth; ++col) {
640 }
641 auto srcPosList = localDestToParentSrc->applyForward(destPosList);
642
643 typename DestImageT::x_iterator destXIter = destImage.row_begin(row);
644 for (int col = 0; col < destWidth; ++col, ++destXIter) {
645 // column index = column + 1 because the first entry in srcPosList is for column -1
646 auto srcPos = srcPosList[col + 1];
647 double relativeArea =
648 computeRelativeArea(srcPos, prevSrcPosList[col], prevSrcPosList[col + 1]);
649
650 if (warpAtOnePoint(destXIter, srcPos, relativeArea,
652 ++numGoodPixels;
653 }
654 } // for col
655 // move points from srcPosList to prevSrcPosList (we don't care about what ends up in srcPosList
656 // because it will be reallocated anyway)
657 swap(srcPosList, prevSrcPosList);
658 } // for row
659 } // if interp
660
661 return numGoodPixels;
662}
663
664template <typename DestImageT, typename SrcImageT>
665int warpCenteredImage(DestImageT &destImage, SrcImageT const &srcImage,
666 lsst::geom::LinearTransform const &linearTransform,
667 lsst::geom::Point2D const &centerPosition, WarpingControl const &control,
668 typename DestImageT::SinglePixel padValue) {
669 // force src and dest to be the same size and xy0
670 if ((destImage.getWidth() != srcImage.getWidth()) || (destImage.getHeight() != srcImage.getHeight()) ||
671 (destImage.getXY0() != srcImage.getXY0())) {
672 std::ostringstream errStream;
673 errStream << "src and dest images must have same size and xy0.";
675 }
676
677 // set the xy0 coords to 0,0 to make life easier
678 SrcImageT srcImageCopy(srcImage, true);
679 srcImageCopy.setXY0(0, 0);
680 destImage.setXY0(0, 0);
681 lsst::geom::Extent2D cLocal =
682 lsst::geom::Extent2D(centerPosition) - lsst::geom::Extent2D(srcImage.getXY0());
683
684 // for the affine transform, the centerPosition will not only get sheared, but also
685 // moved slightly. So we'll include a translation to move it back by an amount
686 // centerPosition - translatedCenterPosition
687 lsst::geom::AffineTransform affTran(linearTransform, cLocal - linearTransform(cLocal));
688 std::shared_ptr<geom::TransformPoint2ToPoint2> affineTransform22 = geom::makeTransform(affTran);
689
690// now warp
691#if 0
692 static float t = 0.0;
693 float t_before = 1.0*clock()/CLOCKS_PER_SEC;
694 int n = warpImage(destImage, srcImageCopy, affTran, control, padValue);
695 float t_after = 1.0*clock()/CLOCKS_PER_SEC;
696 float dt = t_after - t_before;
697 t += dt;
698 std::cout <<srcImage.getWidth()<<"x"<<srcImage.getHeight()<<": "<< dt <<" "<< t <<std::endl;
699#else
700 int n = warpImage(destImage, srcImageCopy, *affineTransform22, control, padValue);
701#endif
702
703 // fix the origin and we're done.
704 destImage.setXY0(srcImage.getXY0());
705
706 return n;
707}
708
709//
710// Explicit instantiations
711//
713// may need to omit default params for EXPOSURE -- original code did that and it worked
714#define EXPOSURE(PIXTYPE) image::Exposure<PIXTYPE, image::MaskPixel, image::VariancePixel>
715#define MASKEDIMAGE(PIXTYPE) image::MaskedImage<PIXTYPE, image::MaskPixel, image::VariancePixel>
716#define IMAGE(PIXTYPE) image::Image<PIXTYPE>
717#define NL /* */
718
719#define INSTANTIATE(DESTIMAGEPIXELT, SRCIMAGEPIXELT) \
720 template int warpCenteredImage( \
721 IMAGE(DESTIMAGEPIXELT) & destImage, IMAGE(SRCIMAGEPIXELT) const &srcImage, \
722 lsst::geom::LinearTransform const &linearTransform, lsst::geom::Point2D const &centerPosition, \
723 WarpingControl const &control, IMAGE(DESTIMAGEPIXELT)::SinglePixel padValue); \
724 NL template int warpCenteredImage( \
725 MASKEDIMAGE(DESTIMAGEPIXELT) & destImage, MASKEDIMAGE(SRCIMAGEPIXELT) const &srcImage, \
726 lsst::geom::LinearTransform const &linearTransform, lsst::geom::Point2D const &centerPosition, \
727 WarpingControl const &control, MASKEDIMAGE(DESTIMAGEPIXELT)::SinglePixel padValue); \
728 NL template int warpImage(IMAGE(DESTIMAGEPIXELT) & destImage, IMAGE(SRCIMAGEPIXELT) const &srcImage, \
729 geom::TransformPoint2ToPoint2 const &srcToDest, WarpingControl const &control, \
730 IMAGE(DESTIMAGEPIXELT)::SinglePixel padValue); \
731 NL template int warpImage(MASKEDIMAGE(DESTIMAGEPIXELT) & destImage, \
732 MASKEDIMAGE(SRCIMAGEPIXELT) const &srcImage, \
733 geom::TransformPoint2ToPoint2 const &srcToDest, WarpingControl const &control, \
734 MASKEDIMAGE(DESTIMAGEPIXELT)::SinglePixel padValue); \
735 NL template int warpImage(IMAGE(DESTIMAGEPIXELT) & destImage, geom::SkyWcs const &destWcs, \
736 IMAGE(SRCIMAGEPIXELT) const &srcImage, geom::SkyWcs const &srcWcs, \
737 WarpingControl const &control, IMAGE(DESTIMAGEPIXELT)::SinglePixel padValue); \
738 NL template int warpImage(MASKEDIMAGE(DESTIMAGEPIXELT) & destImage, geom::SkyWcs const &destWcs, \
739 MASKEDIMAGE(SRCIMAGEPIXELT) const &srcImage, geom::SkyWcs const &srcWcs, \
740 WarpingControl const &control, \
741 MASKEDIMAGE(DESTIMAGEPIXELT)::SinglePixel padValue); \
742 NL template int warpExposure(EXPOSURE(DESTIMAGEPIXELT) & destExposure, \
743 EXPOSURE(SRCIMAGEPIXELT) const &srcExposure, WarpingControl const &control, \
744 EXPOSURE(DESTIMAGEPIXELT)::MaskedImageT::SinglePixel padValue);
745
746INSTANTIATE(double, double)
747INSTANTIATE(double, float)
748INSTANTIATE(double, int)
750INSTANTIATE(float, float)
751INSTANTIATE(float, int)
753INSTANTIATE(int, int)
756} // namespace math
757
758template std::shared_ptr<math::LanczosWarpingKernel> table::io::PersistableFacade<
759 math::LanczosWarpingKernel>::dynamicCast(std::shared_ptr<table::io::Persistable> const &);
760template std::shared_ptr<math::BilinearWarpingKernel> table::io::PersistableFacade<
761 math::BilinearWarpingKernel>::dynamicCast(std::shared_ptr<table::io::Persistable> const &);
762template std::shared_ptr<math::NearestWarpingKernel> table::io::PersistableFacade<
763 math::NearestWarpingKernel>::dynamicCast(std::shared_ptr<table::io::Persistable> const &);
764template std::shared_ptr<math::WarpingControl> table::io::PersistableFacade<
765 math::WarpingControl>::dynamicCast(std::shared_ptr<table::io::Persistable> const &);
766
767} // namespace afw
768} // namespace lsst
int end
#define INSTANTIATE(FROMSYS, TOSYS)
Definition Detector.cc:509
#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
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
table::Schema schema
Definition python.h:134
T back(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.
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
int getWidth() const
Return the Kernel's width.
Definition Kernel.h:224
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.
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.
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.
void setKernelParameter(unsigned int ind, double value) const override
Set one kernel parameter.
SeparableKernel()
Construct an empty spatially invariant SeparableKernel of size 0x0.
void setWarpingKernel(SeparableKernel const &warpingKernel)
set the warping kernel
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.
PersistableFactory(std::string const &name)
Constructor for the factory.
io::OutputArchiveHandle OutputArchiveHandle
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 clock(T... args)
T emplace_back(T... args)
T empty(T... args)
T endl(T... args)
T fabs(T... args)
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:693
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.
CatalogT< BaseRecord > BaseCatalog
Definition fwd.h:72
Extent< double, 2 > Extent2D
Definition Extent.h:400
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:355
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
std::shared_ptr< table::io::Persistable > read(table::io::InputArchive const &archive, table::io::CatalogVector const &catalogs) const override
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

Variable Documentation

◆ cacheSize

table::Key<int> cacheSize

Definition at line 306 of file warpExposure.cc.

◆ growFullMask

table::Key<image::MaskPixel> growFullMask

Definition at line 308 of file warpExposure.cc.

◆ hasMaskKernel

table::Key<table::Flag> hasMaskKernel

Definition at line 304 of file warpExposure.cc.

◆ interpLength

table::Key<int> interpLength

Definition at line 307 of file warpExposure.cc.

◆ maskKernelIndex

table::Key<int> maskKernelIndex

Definition at line 305 of file warpExposure.cc.

◆ order

table::Key<int> order

Definition at line 165 of file warpExposure.cc.

◆ schema

table::Schema schema

Definition at line 164 of file warpExposure.cc.

◆ warpingKernelIndex

table::Key<int> warpingKernelIndex

Definition at line 303 of file warpExposure.cc.