Loading [MathJax]/extensions/tex2jax.js
LSST Applications g04a91732dc+cc8870d3f5,g07dc498a13+5aa0b8792f,g0fba68d861+488cddfaa9,g1409bbee79+5aa0b8792f,g1a7e361dbc+5aa0b8792f,g1fd858c14a+f64bc332a9,g35bb328faa+fcb1d3bbc8,g4d2262a081+b1c1982739,g4d39ba7253+9633a327c1,g4e0f332c67+5d362be553,g53246c7159+fcb1d3bbc8,g60b5630c4e+9633a327c1,g668ecb457e+25d63fd678,g78460c75b0+2f9a1b4bcd,g786e29fd12+cf7ec2a62a,g7b71ed6315+fcb1d3bbc8,g8852436030+8b64ca622a,g89139ef638+5aa0b8792f,g89e1512fd8+37f975783e,g8d6b6b353c+9633a327c1,g9125e01d80+fcb1d3bbc8,g989de1cb63+5aa0b8792f,g9f33ca652e+b196626af7,ga9baa6287d+9633a327c1,gaaedd4e678+5aa0b8792f,gabe3b4be73+1e0a283bba,gb1101e3267+71e32094df,gb58c049af0+f03b321e39,gb90eeb9370+2807b1ad02,gc741bbaa4f+1ae86710ed,gcf25f946ba+8b64ca622a,gd315a588df+a39986a76f,gd6cbbdb0b4+c8606af20c,gd9a9a58781+fcb1d3bbc8,gde0f65d7ad+94dfc458f4,ge278dab8ac+932305ba37,ge82c20c137+76d20ab76d,gfe73954cf8+a1301e4c20,w.2025.11
LSST Data Management Base Package
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
WarpedPsf.cc File Reference

Go to the source code of this file.

Namespaces

namespace  lsst
 
namespace  lsst::afw
 
namespace  lsst::afw::table
 
namespace  lsst::afw::table::io
 
namespace  lsst::meas
 
namespace  lsst::meas::algorithms
 

Functions

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

Function Documentation

◆ read()

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

Definition at line 1 of file WarpedPsf.cc.

22 ://www.lsstcorp.org/LegalNotices/>.
23 */
24
26#include "lsst/geom/Box.h"
27#include "lsst/pex/exceptions.h"
35
36namespace lsst {
37namespace afw {
38namespace table {
39namespace io {
40
41template std::shared_ptr<meas::algorithms::WarpedPsf>
42PersistableFacade<meas::algorithms::WarpedPsf>::dynamicCast(std::shared_ptr<Persistable> const &);
43
44} // namespace io
45} // namespace table
46} // namespace afw
47
48namespace meas {
49namespace algorithms {
50
51namespace {
52
53// TODO: make this routine externally callable and more generic using templates
54// (also useful in e.g. math/offsetImage.cc)
55std::shared_ptr<afw::detection::Psf::Image> zeroPadImage(afw::detection::Psf::Image const &im, int xPad,
56 int yPad) {
57 int nx = im.getWidth();
58 int ny = im.getHeight();
59
60 auto out = std::make_shared<afw::detection::Psf::Image>(nx + 2 * xPad, ny + 2 * yPad);
61 out->setXY0(im.getX0() - xPad, im.getY0() - yPad);
62
63 geom::Box2I box(geom::Point2I(xPad, yPad), geom::Extent2I(nx, ny));
64 out->assign(im, box, afw::image::LOCAL);
65
66 return out;
67}
68
69geom::Box2I computeBBoxFromTransform(geom::Box2I const bbox, geom::AffineTransform const &t) {
70 // This is the maximum scaling I can imagine we'd want - it's approximately what you'd
71 // get from trying to coadd 2"-pixel data (e.g. 2MASS) onto a 0.01"-pixel grid (e.g.
72 // from JWST). Anything beyond that is probably a bug elsewhere in the code, and
73 // catching that can save us from segfaults and catastrophic memory consumption.
74 static const double maxTransformCoeff = 200.0;
75
76 if (t.getLinear().getMatrix().lpNorm<Eigen::Infinity>() > maxTransformCoeff) {
77 throw LSST_EXCEPT(pex::exceptions::RangeError, "Unexpectedly large transform passed to WarpedPsf");
78 }
79
80 // Floating point version of input bbox (expanded to include entire pixels).
81 geom::Box2D in_box_fp(bbox);
82 // Floating point version of output bbox (to be populated as bbox of
83 // transformed points.
84 geom::Box2D out_box_fp;
85 auto const in_corners = in_box_fp.getCorners();
86 for (auto const & in_corner : in_corners) {
87 auto out_corner = t(in_corner);
88 out_box_fp.include(out_corner);
89 }
90
91 // We want to guarantee that the output bbox has odd dimensions, so instead
92 // of using the Box2I converting constructor directly, we start with (0, 0)
93 // and dilate by the floating-point box's half-dimensions.
94 geom::Extent2I out_half_dims = geom::floor(0.5*out_box_fp.getDimensions());
95 geom::Box2I out_box;
96 geom::Point2I out_center(0, 0);
97 out_box.include(out_center);
98 return out_box.dilatedBy(out_half_dims);
99}
100
114std::shared_ptr<afw::detection::Psf::Image> warpAffine(afw::detection::Psf::Image const &im,
115 geom::AffineTransform const &srcToDest,
116 afw::math::WarpingControl const &wc) {
117 std::shared_ptr<afw::geom::TransformPoint2ToPoint2> srcToDestTransform =
118 afw::geom::makeTransform(srcToDest);
119
120 afw::math::SeparableKernel const &kernel = *wc.getWarpingKernel();
121 geom::Point2I const &center = kernel.getCtr();
122 int const xPad = std::max(center.getX(), kernel.getWidth() - center.getX());
123 int const yPad = std::max(center.getY(), kernel.getHeight() - center.getY());
124
125 // allocate output image
126 geom::Box2I bbox = computeBBoxFromTransform(im.getBBox(), srcToDest);
127
129
130 // zero-pad input image
131 std::shared_ptr<afw::detection::Psf::Image> im_padded = zeroPadImage(im, xPad, yPad);
132
133 // warp it!
134 afw::math::warpImage(*ret, *im_padded, *srcToDestTransform, wc, 0.0);
135 return ret;
136}
137
138// A helper struct for logic and intermediate results used by both
139// doComputeKernelImage and doComputeBBox.
140//
141// This linearizes the transform from destination to source coordinates, then
142// computes both a source position that corresponds to a rounded version of the
143// destination position and a transform from source to destination that also
144// undoes the rounding.
145struct PreparedTransforms {
146
147 static PreparedTransforms compute(
148 afw::detection::Psf const & src_psf,
149 afw::geom::TransformPoint2ToPoint2 const & distortion,
150 geom::PointD const & dest_position
151 ) {
152 // Linearize the transform from destination coordinates to source
153 // coordinates to avoid expensive WCS calls.
155 *distortion.inverted(),
156 dest_position
157 );
158 // Instead of calling computeKernelImage on the source PSFs, we want to
159 // call computeImage with the source coordinate version of an
160 // integer-valued destination coordinate; that gives the lower-level PSF
161 // model the responsibility of doing the shifting from integer to
162 // fractional source coordinates, and it may be able to do a better job of
163 // that than we could do here (by e.g. using an internal analytic or
164 // oversampled model). So here's that integer-valued destination position
165 // and the corresponding source position.
166 geom::Point2D rounded_dest_position(
167 std::round(dest_position.getX()),
168 std::round(dest_position.getY())
169 );
170 auto src_position = dest_to_src(rounded_dest_position);
171 // Now we want to warp by the inverse of src_to_dest, but we want the
172 // image to be centered at (0, 0), not rounded_dest_position, so we
173 // compose that trivial shift with the dest_to_src transform.
174 auto src_to_dest = geom::AffineTransform(dest_to_src.inverted());
175 auto unround = geom::AffineTransform(-geom::Extent2D(rounded_dest_position));
176 auto src_to_dest_unrounded = unround * src_to_dest;
177 return PreparedTransforms{src_position, src_to_dest_unrounded};
178 }
179
180 geom::Point2D src_position;
181 geom::AffineTransform src_to_dest_unrounded;
182};
183
184} // namespace
185
186WarpedPsf::WarpedPsf(std::shared_ptr<afw::detection::Psf const> undistortedPsf,
187 std::shared_ptr<afw::geom::TransformPoint2ToPoint2 const> distortion,
188 std::shared_ptr<afw::math::WarpingControl const> control)
189 : ImagePsf(false),
190 _undistortedPsf(undistortedPsf),
191 _distortion(distortion),
192 _warpingControl(control) {
193 _init();
194}
195
196WarpedPsf::WarpedPsf(std::shared_ptr<afw::detection::Psf const> undistortedPsf,
198 std::string const &kernelName, unsigned int cache)
199 : ImagePsf(false),
200 _undistortedPsf(undistortedPsf),
201 _distortion(distortion),
202 _warpingControl(new afw::math::WarpingControl(kernelName, "", cache)) {
203 _init();
204}
205
206void WarpedPsf::_init() {
207 if (!_undistortedPsf) {
208 throw LSST_EXCEPT(pex::exceptions::LogicError,
209 "Undistorted Psf passed to WarpedPsf must not be None/NULL");
210 }
211 if (!_distortion) {
212 throw LSST_EXCEPT(pex::exceptions::LogicError, "Transform passed to WarpedPsf must not be None/NULL");
213 }
214 if (!_warpingControl) {
215 throw LSST_EXCEPT(pex::exceptions::LogicError,
216 "WarpingControl passed to WarpedPsf must not be None/NULL");
217 }
218}
219
220geom::Point2D WarpedPsf::getAveragePosition() const {
221 return _distortion->applyForward(_undistortedPsf->getAveragePosition());
222}
223
224std::shared_ptr<afw::detection::Psf> WarpedPsf::clone() const {
225 return std::make_shared<WarpedPsf>(_undistortedPsf->clone(), _distortion, _warpingControl);
226}
227
228std::shared_ptr<afw::detection::Psf> WarpedPsf::resized(int width, int height) const {
229 // For a given set of requested dimensions and distortion, it is not guaranteed that a
230 // _undistortedPsf would exist to manifest those dimensions after distortion
231 // Not possible to implement with member data currently in WarpedPsf
232 throw LSST_EXCEPT(pex::exceptions::LogicError, "Not Implemented");
233}
234
235std::shared_ptr<afw::detection::Psf::Image> WarpedPsf::doComputeKernelImage(
236 geom::Point2D const &position, afw::image::Color const &color) const {
237 auto prepared_transforms = PreparedTransforms::compute(*_undistortedPsf, *_distortion, position);
238 std::shared_ptr<Image> src_im = _undistortedPsf->computeImage(prepared_transforms.src_position, color);
240 *src_im,
241 prepared_transforms.src_to_dest_unrounded,
242 *_warpingControl
243 );
244
245 double normFactor = 1.0;
246 //
247 // Normalize the output image to sum 1
248 // FIXME defining a member function Image::getSum() would be convenient here and in other places
249 //
250 normFactor = 0.0;
251 for (int y = 0; y != ret->getHeight(); ++y) {
252 afw::detection::Psf::Image::x_iterator imEnd = ret->row_end(y);
253 for (afw::detection::Psf::Image::x_iterator imPtr = ret->row_begin(y); imPtr != imEnd; imPtr++) {
254 normFactor += *imPtr;
255 }
256 }
257 if (normFactor == 0.0) {
258 throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, "psf image has sum 0");
259 }
260 *ret /= normFactor;
261 return ret;
262}
263
264geom::Box2I WarpedPsf::doComputeBBox(geom::Point2D const &position, afw::image::Color const &color) const {
265 auto prepared_transforms = PreparedTransforms::compute(*_undistortedPsf, *_distortion, position);
266 auto src_bbox = _undistortedPsf->computeImageBBox(prepared_transforms.src_position, color);
267 return computeBBoxFromTransform(src_bbox, prepared_transforms.src_to_dest_unrounded);
268}
269
270namespace {
271
272struct PersistenceHelper {
273 afw::table::Schema schema;
274 afw::table::Key<int> psfIndex;
275 afw::table::Key<int> transformIndex;
276 afw::table::Key<int> controlIndex;
277
278 static PersistenceHelper const &get() {
279 static PersistenceHelper const instance;
280 return instance;
281 }
282
283private:
284 PersistenceHelper()
285 : schema(),
286 psfIndex(schema.addField<int>("psfIndex", "archive ID of nested Psf object")),
287 transformIndex(schema.addField<int>("transformIndex", "archive ID of nested Transform object")),
288 controlIndex(
289 schema.addField<int>("controlIndex", "archive ID of nested WarpingControl object")) {}
290};
291
292std::string _getPersistenceName() { return "WarpedPsf"; }
293
294class : public afw::table::io::PersistableFactory {
295public:
297 afw::table::io::InputArchive const &archive,
298 afw::table::io::CatalogVector const &catalogs) const {
299 static PersistenceHelper const &keys = PersistenceHelper::get();
300 LSST_ARCHIVE_ASSERT(catalogs.size() == 1u);
301 LSST_ARCHIVE_ASSERT(catalogs.front().size() == 1u);
302 afw::table::BaseRecord const &record = catalogs.front().front();
303 LSST_ARCHIVE_ASSERT(record.getSchema() == keys.schema);
305 archive.get<afw::detection::Psf>(record.get(keys.psfIndex)),
306 archive.get<afw::geom::TransformPoint2ToPoint2>(record.get(keys.transformIndex)),
307 archive.get<afw::math::WarpingControl>(record.get(keys.controlIndex)));
308 }
309
310 using afw::table::io::PersistableFactory::PersistableFactory;
311} warpedPsfFactory(_getPersistenceName());
312
313} // namespace
314
315std::string WarpedPsf::getPersistenceName() const { return _getPersistenceName(); }
316std::string WarpedPsf::getPythonModule() const { return "lsst.meas.algorithms"; }
317
318void WarpedPsf::write(OutputArchiveHandle &handle) const {
319 static PersistenceHelper const &keys = PersistenceHelper::get();
320 afw::table::BaseCatalog catalog = handle.makeCatalog(keys.schema);
321 std::shared_ptr<afw::table::BaseRecord> record = catalog.addNew();
322
323 record->set(keys.psfIndex, handle.put(_undistortedPsf));
324 record->set(keys.transformIndex, handle.put(_distortion));
325 record->set(keys.controlIndex, handle.put(_warpingControl));
326
327 handle.saveCatalog(catalog);
328}
329
330} // namespace algorithms
331} // namespace meas
332} // namespace lsst
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
Definition Persistable.h:48
A polymorphic base class for representing an image's Point Spread Function.
Definition Psf.h:82
image::Image< Pixel > Image
Image type returned by computeImage.
Definition Psf.h:89
A kernel described by a pair of functions: func(x, y) = colFunc(x) * rowFunc(y)
Definition Kernel.h:860
Parameters to control convolution.
static std::shared_ptr< T > dynamicCast(std::shared_ptr< Persistable > const &ptr)
Dynamically cast a shared_ptr.
An affine coordinate transformation consisting of a linear transformation and an offset.
AffineTransform const inverted() const
Return the inverse transform.
LinearTransform const & getLinear() const noexcept
A floating-point coordinate rectangle geometry.
Definition Box.h:413
Extent2D const getDimensions() const noexcept
1-d interval accessors
Definition Box.h:528
void include(Point2D const &point) noexcept
Expand this to ensure that this->contains(point).
Definition Box.cc:380
std::vector< Point2D > getCorners() const
Get the corner points.
Definition Box.cc:496
An integer coordinate rectangle.
Definition Box.h:55
void include(Point2I const &point)
Expand this to ensure that this->contains(point).
Definition Box.cc:152
Box2I dilatedBy(Extent const &buffer) const
Increase the size of the box by the given amount(s) on all sides (returning a new object).
Definition Box.cc:213
Matrix const & getMatrix() const noexcept
Reports when the result of an operation cannot be represented by the destination type.
Definition Runtime.h:115
T make_shared(T... args)
T max(T... args)
std::shared_ptr< TransformPoint2ToPoint2 > makeTransform(lsst::geom::AffineTransform const &affine)
Wrap an lsst::geom::AffineTransform as a Transform.
Transform< Point2Endpoint, Point2Endpoint > TransformPoint2ToPoint2
Definition Transform.h:300
lsst::geom::AffineTransform linearizeTransform(TransformPoint2ToPoint2 const &original, lsst::geom::Point2D const &inPoint)
Approximate a Transform by its local linearization.
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.
Extent< int, N > floor(Extent< double, N > const &input) noexcept
Return the component-wise floor (round towards more negative).
Definition Extent.cc:109
Point< double, 2 > PointD
Definition Point.h:323
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 round(T... args)
virtual std::shared_ptr< afw::table::io::Persistable > read(afw::table::io::InputArchive const &archive, afw::table::io::CatalogVector const &catalogs) const
Definition WarpedPsf.cc:1