LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
WarpedPsf.cc
Go to the documentation of this file.
1// -*- lsst-c++ -*-
2
3/*
4 * LSST Data Management System
5 * Copyright 2008, 2009, 2010 LSST Corporation.
6 *
7 * This product includes software developed by the
8 * LSST Project (http://www.lsst.org/).
9 *
10 * This program is free software: you can redistribute it and/or modify
11 * it under the terms of the GNU General Public License as published by
12 * the Free Software Foundation, either version 3 of the License, or
13 * (at your option) any later version.
14 *
15 * This program is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 * GNU General Public License for more details.
19 *
20 * You should have received a copy of the LSST License Statement and
21 * the GNU General Public License along with this program. If not,
22 * see <http://www.lsstcorp.org/LegalNotices/>.
23 */
24
26#include "lsst/geom/Box.h"
27#include "lsst/pex/exceptions.h"
35
36namespace lsst {
37namespace afw {
38namespace table {
39namespace io {
40
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)
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
115 geom::AffineTransform const &srcToDest,
116 afw::math::WarpingControl const &wc) {
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
128 auto ret = std::make_shared<afw::detection::Psf::Image>(bbox);
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
182};
183
184} // namespace
185
186WarpedPsf::WarpedPsf(std::shared_ptr<afw::detection::Psf const> undistortedPsf,
189 : ImagePsf(false),
190 _undistortedPsf(undistortedPsf),
191 _distortion(distortion),
192 _warpingControl(control) {
193 _init();
194}
195
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) {
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) {
216 "WarpingControl passed to WarpedPsf must not be None/NULL");
217 }
218}
219
221 return _distortion->applyForward(_undistortedPsf->getAveragePosition());
222}
223
225 return std::make_shared<WarpedPsf>(_undistortedPsf->clone(), _distortion, _warpingControl);
226}
227
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
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
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")),
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);
304 return std::make_shared<WarpedPsf>(
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
311} warpedPsfFactory(_getPersistenceName());
312
313} // namespace
314
315std::string WarpedPsf::getPersistenceName() const { return _getPersistenceName(); }
316std::string WarpedPsf::getPythonModule() const { return "lsst.meas.algorithms"; }
317
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
AmpInfoBoxKey bbox
Definition Amplifier.cc:117
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
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
image::Image< Pixel > Image
Image type returned by computeImage.
Definition Psf.h:83
Describe the colour of a source.
Definition Color.h:25
typename _view_t::x_iterator x_iterator
An iterator for traversing the pixels in a row.
Definition ImageBase.h:133
An object passed to Persistable::write to allow it to persist itself.
void saveCatalog(BaseCatalog const &catalog)
Save a catalog in the archive.
BaseCatalog makeCatalog(Schema const &schema)
Return a new, empty catalog with the given schema.
int put(Persistable const *obj, bool permissive=false)
Save an object to the archive and return a unique ID that can be used to retrieve it from an InputArc...
static std::shared_ptr< T > dynamicCast(std::shared_ptr< Persistable > const &ptr)
Dynamically cast a shared_ptr.
PersistableFactory(std::string const &name)
Constructor for the factory.
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
An intermediate base class for Psfs that use an image representation.
Definition ImagePsf.h:40
void write(OutputArchiveHandle &handle) const override
Write the object to one or more catalogs.
Definition WarpedPsf.cc:318
std::shared_ptr< afw::detection::Psf > resized(int width, int height) const override
Return a clone with specified kernel dimensions.
Definition WarpedPsf.cc:228
geom::Point2D getAveragePosition() const override
Return the average of the positions of the stars that went into this Psf.
Definition WarpedPsf.cc:220
geom::Box2I doComputeBBox(geom::Point2D const &position, afw::image::Color const &color) const override
Definition WarpedPsf.cc:264
std::shared_ptr< afw::detection::Psf::Image > doComputeKernelImage(geom::Point2D const &position, afw::image::Color const &color) const override
These virtual member functions are private, not protected, because we only want derived classes to im...
Definition WarpedPsf.cc:235
std::shared_ptr< afw::geom::TransformPoint2ToPoint2 const > _distortion
Definition WarpedPsf.h:94
std::shared_ptr< afw::detection::Psf const > _undistortedPsf
Definition WarpedPsf.h:93
WarpedPsf(std::shared_ptr< afw::detection::Psf const > undistortedPsf, std::shared_ptr< afw::geom::TransformPoint2ToPoint2 const > distortion, std::shared_ptr< afw::math::WarpingControl const > control)
Construct WarpedPsf from unwarped psf and distortion.
Definition WarpedPsf.cc:186
std::string getPersistenceName() const override
Return the unique name used to persist this object and look up its factory.
Definition WarpedPsf.cc:315
std::shared_ptr< afw::detection::Psf > clone() const override
Polymorphic deep copy. Usually unnecessary, as Psfs are immutable.
Definition WarpedPsf.cc:224
std::string getPythonModule() const override
Return the fully-qualified Python module that should be imported to guarantee that its factory is reg...
Definition WarpedPsf.cc:316
Reports invalid arguments.
Definition Runtime.h:66
Reports errors in the logical structure of the program.
Definition Runtime.h:46
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
T round(T... args)
afw::table::Key< int > controlIndex
Definition WarpedPsf.cc:276
afw::table::Key< int > transformIndex
Definition WarpedPsf.cc:275
afw::table::Key< int > psfIndex
Definition WarpedPsf.cc:274
geom::Point2D src_position
Definition WarpedPsf.cc:180
geom::AffineTransform src_to_dest_unrounded
Definition WarpedPsf.cc:181
std::shared_ptr< table::io::Persistable > read(table::io::InputArchive const &archive, table::io::CatalogVector const &catalogs) const override