LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
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"
30 #include "lsst/afw/image/Image.h"
35 
36 namespace lsst {
37 namespace afw {
38 namespace table {
39 namespace io {
40 
43 
44 } // namespace io
45 } // namespace table
46 } // namespace afw
47 
48 namespace meas {
49 namespace algorithms {
50 
51 namespace {
52 
53 inline double min4(double a, double b, double c, double d) {
54  return std::min(std::min(a, b), std::min(c, d));
55 }
56 
57 inline double max4(double a, double b, double c, double d) {
58  return std::max(std::max(a, b), std::max(c, d));
59 }
60 
61 // TODO: make this routine externally callable and more generic using templates
62 // (also useful in e.g. math/offsetImage.cc)
64  int yPad) {
65  int nx = im.getWidth();
66  int ny = im.getHeight();
67 
68  auto out = std::make_shared<afw::detection::Psf::Image>(nx + 2 * xPad, ny + 2 * yPad);
69  out->setXY0(im.getX0() - xPad, im.getY0() - yPad);
70 
71  geom::Box2I box(geom::Point2I(xPad, yPad), geom::Extent2I(nx, ny));
72  out->assign(im, box, afw::image::LOCAL);
73 
74  return out;
75 }
76 
77 geom::Box2I computeBBoxFromTransform(geom::Box2I const bbox, geom::AffineTransform const &t) {
78  static const int dst_padding = 0;
79 
80  // This is the maximum scaling I can imagine we'd want - it's approximately what you'd
81  // get from trying to coadd 2"-pixel data (e.g. 2MASS) onto a 0.01"-pixel grid (e.g.
82  // from JWST). Anything beyond that is probably a bug elsewhere in the code, and
83  // catching that can save us from segfaults and catastrophic memory consumption.
84  static const double maxTransformCoeff = 200.0;
85 
86  if (t.getLinear().getMatrix().lpNorm<Eigen::Infinity>() > maxTransformCoeff) {
87  throw LSST_EXCEPT(pex::exceptions::RangeError, "Unexpectedly large transform passed to WarpedPsf");
88  }
89 
90  // min/max coordinate values in input image
91  int in_xlo = bbox.getMinX();
92  int in_xhi = bbox.getMinX() + bbox.getWidth() - 1;
93  int in_ylo = bbox.getMinY();
94  int in_yhi = bbox.getMinY() + bbox.getHeight() - 1;
95 
96  // corners of output image
97  geom::Point2D c00 = t(geom::Point2D(in_xlo, in_ylo));
98  geom::Point2D c01 = t(geom::Point2D(in_xlo, in_yhi));
99  geom::Point2D c10 = t(geom::Point2D(in_xhi, in_ylo));
100  geom::Point2D c11 = t(geom::Point2D(in_xhi, in_yhi));
101 
102  //
103  // bounding box for output image
104  //
105  int out_xlo = floor(min4(c00.getX(), c01.getX(), c10.getX(), c11.getX())) - dst_padding;
106  int out_ylo = floor(min4(c00.getY(), c01.getY(), c10.getY(), c11.getY())) - dst_padding;
107  int out_xhi = ceil(max4(c00.getX(), c01.getX(), c10.getX(), c11.getX())) + dst_padding;
108  int out_yhi = ceil(max4(c00.getY(), c01.getY(), c10.getY(), c11.getY())) + dst_padding;
109 
110  geom::Box2I ret = geom::Box2I(geom::Point2I(out_xlo, out_ylo),
111  geom::Extent2I(out_xhi - out_xlo + 1, out_yhi - out_ylo + 1));
112  return ret;
113 }
114 
129  geom::AffineTransform const &srcToDest,
130  afw::math::WarpingControl const &wc) {
132  afw::geom::makeTransform(srcToDest);
133 
134  afw::math::SeparableKernel const &kernel = *wc.getWarpingKernel();
135  geom::Point2I const &center = kernel.getCtr();
136  int const xPad = std::max(center.getX(), kernel.getWidth() - center.getX());
137  int const yPad = std::max(center.getY(), kernel.getHeight() - center.getY());
138 
139  // allocate output image
140  geom::Box2I bbox = computeBBoxFromTransform(im.getBBox(), srcToDest);
141  auto ret = std::make_shared<afw::detection::Psf::Image>(bbox);
142 
143  // zero-pad input image
144  std::shared_ptr<afw::detection::Psf::Image> im_padded = zeroPadImage(im, xPad, yPad);
145 
146  // warp it!
147  afw::math::warpImage(*ret, *im_padded, *srcToDestTransform, wc, 0.0);
148  return ret;
149 }
150 
151 } // namespace
152 
153 WarpedPsf::WarpedPsf(std::shared_ptr<afw::detection::Psf const> undistortedPsf,
156  : ImagePsf(false),
157  _undistortedPsf(undistortedPsf),
158  _distortion(distortion),
159  _warpingControl(control) {
160  _init();
161 }
162 
165  std::string const &kernelName, unsigned int cache)
166  : ImagePsf(false),
167  _undistortedPsf(undistortedPsf),
168  _distortion(distortion),
169  _warpingControl(new afw::math::WarpingControl(kernelName, "", cache)) {
170  _init();
171 }
172 
173 void WarpedPsf::_init() {
174  if (!_undistortedPsf) {
176  "Undistorted Psf passed to WarpedPsf must not be None/NULL");
177  }
178  if (!_distortion) {
179  throw LSST_EXCEPT(pex::exceptions::LogicError, "Transform passed to WarpedPsf must not be None/NULL");
180  }
181  if (!_warpingControl) {
183  "WarpingControl passed to WarpedPsf must not be None/NULL");
184  }
185 }
186 
188  return _distortion->applyForward(_undistortedPsf->getAveragePosition());
189 }
190 
192  return std::make_shared<WarpedPsf>(_undistortedPsf->clone(), _distortion, _warpingControl);
193 }
194 
196  // For a given set of requested dimensions and distortion, it is not guaranteed that a
197  // _undistortedPsf would exist to manifest those dimensions after distortion
198  // Not possible to implement with member data currently in WarpedPsf
199  throw LSST_EXCEPT(pex::exceptions::LogicError, "Not Implemented");
200 }
201 
203  geom::Point2D const &position, afw::image::Color const &color) const {
205  geom::Point2D tp = t(position);
206 
207  std::shared_ptr<Image> im = _undistortedPsf->computeKernelImage(tp, color);
208 
209  // Go to the warped coordinate system with 'p' at the origin
210  auto srcToDest = geom::AffineTransform(t.inverted().getLinear());
211  std::shared_ptr<afw::detection::Psf::Psf::Image> ret = warpAffine(*im, srcToDest, *_warpingControl);
212 
213  double normFactor = 1.0;
214  //
215  // Normalize the output image to sum 1
216  // FIXME defining a member function Image::getSum() would be convenient here and in other places
217  //
218  normFactor = 0.0;
219  for (int y = 0; y != ret->getHeight(); ++y) {
220  afw::detection::Psf::Image::x_iterator imEnd = ret->row_end(y);
221  for (afw::detection::Psf::Image::x_iterator imPtr = ret->row_begin(y); imPtr != imEnd; imPtr++) {
222  normFactor += *imPtr;
223  }
224  }
225  if (normFactor == 0.0) {
226  throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, "psf image has sum 0");
227  }
228  *ret /= normFactor;
229  return ret;
230 }
231 
232 geom::Box2I WarpedPsf::doComputeBBox(geom::Point2D const &position, afw::image::Color const &color) const {
234  geom::Point2D tp = t(position);
235  geom::Box2I bboxUndistorted = _undistortedPsf->computeBBox(tp, color);
236  geom::Box2I ret =
237  computeBBoxFromTransform(bboxUndistorted, geom::AffineTransform(t.inverted().getLinear()));
238  return ret;
239 }
240 
241 namespace {
242 
243 struct PersistenceHelper {
244  afw::table::Schema schema;
245  afw::table::Key<int> psfIndex;
246  afw::table::Key<int> transformIndex;
247  afw::table::Key<int> controlIndex;
248 
249  static PersistenceHelper const &get() {
250  static PersistenceHelper const instance;
251  return instance;
252  }
253 
254 private:
255  PersistenceHelper()
256  : schema(),
257  psfIndex(schema.addField<int>("psfIndex", "archive ID of nested Psf object")),
258  transformIndex(schema.addField<int>("transformIndex", "archive ID of nested Transform object")),
259  controlIndex(
260  schema.addField<int>("controlIndex", "archive ID of nested WarpingControl object")) {}
261 };
262 
263 std::string _getPersistenceName() { return "WarpedPsf"; }
264 
265 class : public afw::table::io::PersistableFactory {
266 public:
268  afw::table::io::InputArchive const &archive,
269  afw::table::io::CatalogVector const &catalogs) const {
270  static PersistenceHelper const &keys = PersistenceHelper::get();
271  LSST_ARCHIVE_ASSERT(catalogs.size() == 1u);
272  LSST_ARCHIVE_ASSERT(catalogs.front().size() == 1u);
273  afw::table::BaseRecord const &record = catalogs.front().front();
274  LSST_ARCHIVE_ASSERT(record.getSchema() == keys.schema);
275  return std::make_shared<WarpedPsf>(
276  archive.get<afw::detection::Psf>(record.get(keys.psfIndex)),
277  archive.get<afw::geom::TransformPoint2ToPoint2>(record.get(keys.transformIndex)),
278  archive.get<afw::math::WarpingControl>(record.get(keys.controlIndex)));
279  }
280 
282 } warpedPsfFactory(_getPersistenceName());
283 
284 } // namespace
285 
286 std::string WarpedPsf::getPersistenceName() const { return _getPersistenceName(); }
287 std::string WarpedPsf::getPythonModule() const { return "lsst.meas.algorithms"; }
288 
290  static PersistenceHelper const &keys = PersistenceHelper::get();
291  afw::table::BaseCatalog catalog = handle.makeCatalog(keys.schema);
293 
294  record->set(keys.psfIndex, handle.put(_undistortedPsf));
295  record->set(keys.transformIndex, handle.put(_distortion));
296  record->set(keys.controlIndex, handle.put(_warpingControl));
297 
298  handle.saveCatalog(catalog);
299 }
300 
301 } // namespace algorithms
302 } // namespace meas
303 } // 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
table::Key< int > b
table::Key< int > a
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
Definition: Persistable.h:48
image::Image< Pixel > Image
Image type returned by computeImage.
Definition: Psf.h:83
Describe the colour of a source.
Definition: Color.h:26
typename _view_t::x_iterator x_iterator
An iterator for traversing the pixels in a row.
Definition: ImageBase.h:133
std::shared_ptr< RecordT > addNew()
Create a new record, add it to the end of the catalog, and return a pointer to it.
Definition: Catalog.h:490
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.
Definition: Persistable.cc:18
PersistableFactory(std::string const &name)
Constructor for the factory.
Definition: Persistable.cc:74
An affine coordinate transformation consisting of a linear transformation and an offset.
AffineTransform const inverted() const
Return the inverse transform.
LinearTransform const & getLinear() const noexcept
An integer coordinate rectangle.
Definition: Box.h:55
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:289
std::shared_ptr< afw::detection::Psf > resized(int width, int height) const override
Return a clone with specified kernel dimensions.
Definition: WarpedPsf.cc:195
geom::Point2D getAveragePosition() const override
Return the average of the positions of the stars that went into this Psf.
Definition: WarpedPsf.cc:187
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:202
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:153
std::string getPersistenceName() const override
Return the unique name used to persist this object and look up its factory.
Definition: WarpedPsf.cc:286
std::shared_ptr< afw::detection::Psf > clone() const override
Polymorphic deep copy. Usually unnecessary, as Psfs are immutable.
Definition: WarpedPsf.cc:191
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:287
Reports invalid arguments.
Definition: Runtime.h:66
Reports errors in the logical structure of the program.
Definition: Runtime.h:46
T max(T... args)
T min(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
Extent< int, N > ceil(Extent< double, N > const &input) noexcept
Return the component-wise ceil (round towards more positive).
Definition: Extent.cc:118
A base class for image defects.
afw::table::Key< int > controlIndex
Definition: WarpedPsf.cc:247
afw::table::Key< int > transformIndex
Definition: WarpedPsf.cc:246
afw::table::Key< int > psfIndex
Definition: WarpedPsf.cc:245
afw::table::Schema schema
Definition: WarpedPsf.cc:244
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