Loading [MathJax]/extensions/tex2jax.js
LSST Applications g04a91732dc+a777afbe81,g07dc498a13+7e3c5f68a2,g12483e3c20+0145ec33cd,g1409bbee79+7e3c5f68a2,g1a7e361dbc+7e3c5f68a2,g1fd858c14a+9f35e23ec3,g35bb328faa+fcb1d3bbc8,g3ad4f90e5c+0145ec33cd,g3bd4b5ce2c+cbf1bea503,g4e0f332c67+5d362be553,g53246c7159+fcb1d3bbc8,g5477a8d5ce+db04660fe6,g60b5630c4e+0145ec33cd,g623d845a50+0145ec33cd,g6f0c2978f1+3526b51a37,g75b6c65c88+d54b601591,g78460c75b0+2f9a1b4bcd,g786e29fd12+cf7ec2a62a,g7b71ed6315+fcb1d3bbc8,g8852436030+4639f750a5,g89139ef638+7e3c5f68a2,g9125e01d80+fcb1d3bbc8,g919ac25b3e+6220c5324a,g95236ca021+f7a31438ed,g989de1cb63+7e3c5f68a2,g9f33ca652e+2d6fa11d35,gaaedd4e678+7e3c5f68a2,gabe3b4be73+1e0a283bba,gb1101e3267+4a428ef779,gb4a253aaf5+0122250889,gb58c049af0+f03b321e39,gc99c83e5f0+76d20ab76d,gcf25f946ba+4639f750a5,gd6cbbdb0b4+c8606af20c,gde0f65d7ad+3d8a3b7e46,ge278dab8ac+932305ba37,gf795337580+03b96afe58,gfba249425e+fcb1d3bbc8,w.2025.08
LSST Data Management Base Package
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
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.
 
swap (T... args)
 

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
52namespace pexExcept = lsst::pex::exceptions;
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) {
74 throw LSST_EXCEPT(pexExcept::InvalidParameterError,
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) {
82 throw LSST_EXCEPT(pexExcept::InvalidParameterError,
83 "bad coordinate in WarpingKernel::setKernelParameter()");
84 }
85 } else if (ctr == (size + 1) / 2) {
86 if (value < -1 - 1e-6 || value > 1e-6) {
87 throw LSST_EXCEPT(pexExcept::InvalidParameterError,
88 "bad coordinate in WarpingKernel::setKernelParameter()");
89 }
90 } else {
91 throw LSST_EXCEPT(pexExcept::InvalidParameterError,
92 "bad ctr value in WarpingKernel::setKernelParameter()");
93 }
94}
95
96std::shared_ptr<Kernel> LanczosWarpingKernel::clone() const {
97 return std::shared_ptr<Kernel>(new LanczosWarpingKernel(this->getOrder()));
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
107std::shared_ptr<Kernel> BilinearWarpingKernel::clone() const {
108 return std::shared_ptr<Kernel>(new BilinearWarpingKernel());
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
132std::string BilinearWarpingKernel::BilinearFunction1::toString(std::string const &prefix) const {
133 std::ostringstream os;
134 os << "_BilinearFunction1: ";
135 os << Function1<Kernel::Pixel>::toString(prefix);
136 return os.str();
137}
138
139std::shared_ptr<Kernel> NearestWarpingKernel::clone() const {
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
153std::string NearestWarpingKernel::NearestFunction1::toString(std::string const &prefix) const {
154 std::ostringstream os;
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 {
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
223std::shared_ptr<SeparableKernel> makeWarpingKernel(std::string name) {
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
240std::shared_ptr<SeparableKernel> WarpingControl::getWarpingKernel() const {
241 if (_warpingKernelPtr->getCacheSize() != _cacheSize) {
242 _warpingKernelPtr->computeCache(_cacheSize);
243 }
244 return _warpingKernelPtr;
245};
246
247void WarpingControl::setWarpingKernelName(std::string const &warpingKernelName) {
248 std::shared_ptr<SeparableKernel> warpingKernelPtr(makeWarpingKernel(warpingKernelName));
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
261std::shared_ptr<SeparableKernel> WarpingControl::getMaskWarpingKernel() const {
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(),
322 warpingKernelIndex(
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));
461 } catch (lsst::pex::exceptions::InvalidParameterError const &) {
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 =
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
544 std::vector<lsst::geom::Point2D> endColPosList;
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
579 std::vector<lsst::geom::Point2D> destRowPosList;
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
629 std::vector<lsst::geom::Point2D> destPosList;
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) {
639 destPosList.emplace_back(lsst::geom::Point2D(col, row));
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.";
674 throw LSST_EXCEPT(pexExcept::InvalidParameterError, errStream.str());
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)
749INSTANTIATE(double, std::uint16_t)
750INSTANTIATE(float, float)
751INSTANTIATE(float, int)
752INSTANTIATE(float, std::uint16_t)
753INSTANTIATE(int, int)
754INSTANTIATE(std::uint16_t, std::uint16_t)
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
#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.
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
Definition Persistable.h:48
T back(T... args)
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.
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.
A kernel described by a pair of functions: func(x, y) = colFunc(x) * rowFunc(y)
Definition Kernel.h:860
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
bool contains(Point2I const &point) const noexcept
Return true if the box contains the point.
Definition Box.cc:114
T clear(T... args)
T clock(T... args)
T copy(T... args)
T emplace_back(T... args)
T empty(T... args)
T endl(T... args)
T fabs(T... args)
T make_shared(T... args)
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: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
Point< double, 2 > Point2D
Definition Point.h:324
Extent< int, 2 > Extent2I
Definition Extent.h:397
Point< int, 2 > Point2I
Definition Point.h:321
T static_pointer_cast(T... args)
T push_back(T... args)
T regex_match(T... args)
T reserve(T... args)
T size(T... args)
T stoi(T... args)
T str(T... args)
typename ImageT::image_category image_category
Definition ImageBase.h:67
T swap(T... args)
T swap(T... args)
std::shared_ptr< table::io::Persistable > read(table::io::InputArchive const &archive, table::io::CatalogVector const &catalogs) const override

◆ swap()

T std::swap ( T... args)