LSST Applications  21.0.0+04719a4bac,21.0.0-1-ga51b5d4+f5e6047307,21.0.0-11-g2b59f77+a9c1acf22d,21.0.0-11-ga42c5b2+86977b0b17,21.0.0-12-gf4ce030+76814010d2,21.0.0-13-g1721dae+760e7a6536,21.0.0-13-g3a573fe+768d78a30a,21.0.0-15-g5a7caf0+f21cbc5713,21.0.0-16-g0fb55c1+b60e2d390c,21.0.0-19-g4cded4ca+71a93a33c0,21.0.0-2-g103fe59+bb20972958,21.0.0-2-g45278ab+04719a4bac,21.0.0-2-g5242d73+3ad5d60fb1,21.0.0-2-g7f82c8f+8babb168e8,21.0.0-2-g8f08a60+06509c8b61,21.0.0-2-g8faa9b5+616205b9df,21.0.0-2-ga326454+8babb168e8,21.0.0-2-gde069b7+5e4aea9c2f,21.0.0-2-gecfae73+1d3a86e577,21.0.0-2-gfc62afb+3ad5d60fb1,21.0.0-25-g1d57be3cd+e73869a214,21.0.0-3-g357aad2+ed88757d29,21.0.0-3-g4a4ce7f+3ad5d60fb1,21.0.0-3-g4be5c26+3ad5d60fb1,21.0.0-3-g65f322c+e0b24896a3,21.0.0-3-g7d9da8d+616205b9df,21.0.0-3-ge02ed75+a9c1acf22d,21.0.0-4-g591bb35+a9c1acf22d,21.0.0-4-g65b4814+b60e2d390c,21.0.0-4-gccdca77+0de219a2bc,21.0.0-4-ge8a399c+6c55c39e83,21.0.0-5-gd00fb1e+05fce91b99,21.0.0-6-gc675373+3ad5d60fb1,21.0.0-64-g1122c245+4fb2b8f86e,21.0.0-7-g04766d7+cd19d05db2,21.0.0-7-gdf92d54+04719a4bac,21.0.0-8-g5674e7b+d1bd76f71f,master-gac4afde19b+a9c1acf22d,w.2021.13
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"
31 
32 namespace lsst {
33 namespace meas {
34 namespace algorithms {
35 
36 namespace {
37 
38 inline double min4(double a, double b, double c, double d) {
39  return std::min(std::min(a, b), std::min(c, d));
40 }
41 
42 inline double max4(double a, double b, double c, double d) {
43  return std::max(std::max(a, b), std::max(c, d));
44 }
45 
46 // TODO: make this routine externally callable and more generic using templates
47 // (also useful in e.g. math/offsetImage.cc)
48 PTR(afw::detection::Psf::Image) zeroPadImage(afw::detection::Psf::Image const &im, int xPad, int yPad) {
49  int nx = im.getWidth();
50  int ny = im.getHeight();
51 
53  out = std::make_shared<afw::detection::Psf::Image>(nx + 2 * xPad, ny + 2 * yPad);
54  out->setXY0(im.getX0() - xPad, im.getY0() - yPad);
55 
56  geom::Box2I box(geom::Point2I(xPad, yPad), geom::Extent2I(nx, ny));
57  out->assign(im, box, afw::image::LOCAL);
58 
59  return out;
60 }
61 
62 geom::Box2I computeBBoxFromTransform(geom::Box2I const bbox, geom::AffineTransform const &t) {
63  static const int dst_padding = 0;
64 
65  // This is the maximum scaling I can imagine we'd want - it's approximately what you'd
66  // get from trying to coadd 2"-pixel data (e.g. 2MASS) onto a 0.01"-pixel grid (e.g.
67  // from JWST). Anything beyond that is probably a bug elsewhere in the code, and
68  // catching that can save us from segfaults and catastrophic memory consumption.
69  static const double maxTransformCoeff = 200.0;
70 
71  if (t.getLinear().getMatrix().lpNorm<Eigen::Infinity>() > maxTransformCoeff) {
72  throw LSST_EXCEPT(pex::exceptions::RangeError, "Unexpectedly large transform passed to WarpedPsf");
73  }
74 
75  // min/max coordinate values in input image
76  int in_xlo = bbox.getMinX();
77  int in_xhi = bbox.getMinX() + bbox.getWidth() - 1;
78  int in_ylo = bbox.getMinY();
79  int in_yhi = bbox.getMinY() + bbox.getHeight() - 1;
80 
81  // corners of output image
82  geom::Point2D c00 = t(geom::Point2D(in_xlo, in_ylo));
83  geom::Point2D c01 = t(geom::Point2D(in_xlo, in_yhi));
84  geom::Point2D c10 = t(geom::Point2D(in_xhi, in_ylo));
85  geom::Point2D c11 = t(geom::Point2D(in_xhi, in_yhi));
86 
87  //
88  // bounding box for output image
89  //
90  int out_xlo = floor(min4(c00.getX(), c01.getX(), c10.getX(), c11.getX())) - dst_padding;
91  int out_ylo = floor(min4(c00.getY(), c01.getY(), c10.getY(), c11.getY())) - dst_padding;
92  int out_xhi = ceil(max4(c00.getX(), c01.getX(), c10.getX(), c11.getX())) + dst_padding;
93  int out_yhi = ceil(max4(c00.getY(), c01.getY(), c10.getY(), c11.getY())) + dst_padding;
94 
95  geom::Box2I ret = geom::Box2I(geom::Point2I(out_xlo, out_ylo),
96  geom::Extent2I(out_xhi - out_xlo + 1, out_yhi - out_ylo + 1));
97  return ret;
98 }
99 
114 warpAffine(afw::detection::Psf::Image const &im, geom::AffineTransform const &srcToDest,
115  afw::math::WarpingControl const &wc) {
117  afw::geom::makeTransform(srcToDest);
118 
119  afw::math::SeparableKernel const &kernel = *wc.getWarpingKernel();
120  geom::Point2I const &center = kernel.getCtr();
121  int const xPad = std::max(center.getX(), kernel.getWidth() - center.getX());
122  int const yPad = std::max(center.getY(), kernel.getHeight() - center.getY());
123 
124  // allocate output image
125  geom::Box2I bbox = computeBBoxFromTransform(im.getBBox(), srcToDest);
126  PTR(afw::detection::Psf::Image) ret = std::make_shared<afw::detection::Psf::Image>(bbox);
127 
128  // zero-pad input image
129  PTR(afw::detection::Psf::Image) im_padded = zeroPadImage(im, xPad, yPad);
130 
131  // warp it!
132  afw::math::warpImage(*ret, *im_padded, *srcToDestTransform, wc, 0.0);
133  return ret;
134 }
135 
136 } // namespace
137 
139  PTR(afw::geom::TransformPoint2ToPoint2 const) distortion,
141  : ImagePsf(false),
142  _undistortedPsf(undistortedPsf),
143  _distortion(distortion),
144  _warpingControl(control) {
145  _init();
146 }
147 
149  PTR(afw::geom::TransformPoint2ToPoint2 const) distortion, std::string const &kernelName,
150  unsigned int cache)
151  : ImagePsf(false),
152  _undistortedPsf(undistortedPsf),
153  _distortion(distortion),
154  _warpingControl(new afw::math::WarpingControl(kernelName, "", cache)) {
155  _init();
156 }
157 
158 void WarpedPsf::_init() {
159  if (!_undistortedPsf) {
161  "Undistorted Psf passed to WarpedPsf must not be None/NULL");
162  }
163  if (!_distortion) {
164  throw LSST_EXCEPT(pex::exceptions::LogicError, "Transform passed to WarpedPsf must not be None/NULL");
165  }
166  if (!_warpingControl) {
168  "WarpingControl passed to WarpedPsf must not be None/NULL");
169  }
170 }
171 
173  return _distortion->applyForward(_undistortedPsf->getAveragePosition());
174 }
175 
177  return std::make_shared<WarpedPsf>(_undistortedPsf->clone(), _distortion, _warpingControl);
178 }
179 
180 PTR(afw::detection::Psf) WarpedPsf::resized(int width, int height) const {
181  // For a given set of requested dimensions and distortion, it is not guaranteed that a
182  // _undistortedPsf would exist to manifest those dimensions after distortion
183  // Not possible to implement with member data currently in WarpedPsf
184  throw LSST_EXCEPT(pex::exceptions::LogicError, "Not Implemented");
185 }
186 
188 WarpedPsf::doComputeKernelImage(geom::Point2D const &position, afw::image::Color const &color) const {
189  geom::AffineTransform t = afw::geom::linearizeTransform(*_distortion->inverted(), position);
190  geom::Point2D tp = t(position);
191 
192  PTR(Image) im = _undistortedPsf->computeKernelImage(tp, color);
193 
194  // Go to the warped coordinate system with 'p' at the origin
195  auto srcToDest = geom::AffineTransform(t.inverted().getLinear());
196  PTR(afw::detection::Psf::Psf::Image) ret = warpAffine(*im, srcToDest, *_warpingControl);
197 
198  double normFactor = 1.0;
199  //
200  // Normalize the output image to sum 1
201  // FIXME defining a member function Image::getSum() would be convenient here and in other places
202  //
203  normFactor = 0.0;
204  for (int y = 0; y != ret->getHeight(); ++y) {
205  afw::detection::Psf::Image::x_iterator imEnd = ret->row_end(y);
206  for (afw::detection::Psf::Image::x_iterator imPtr = ret->row_begin(y); imPtr != imEnd; imPtr++) {
207  normFactor += *imPtr;
208  }
209  }
210  if (normFactor == 0.0) {
211  throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, "psf image has sum 0");
212  }
213  *ret /= normFactor;
214  return ret;
215 }
216 
219  geom::Point2D tp = t(position);
220  geom::Box2I bboxUndistorted = _undistortedPsf->computeBBox(tp, color);
221  geom::Box2I ret =
222  computeBBoxFromTransform(bboxUndistorted, geom::AffineTransform(t.inverted().getLinear()));
223  return ret;
224 }
225 
226 } // namespace algorithms
227 } // namespace meas
228 } // 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:49
table::Key< int > b
table::Key< int > a
#define PTR(...)
Definition: base.h:41
#define CONST_PTR(...)
A shared pointer to a const object.
Definition: base.h:47
A polymorphic base class for representing an image's Point Spread Function.
Definition: Psf.h:76
image::Image< Pixel > Image
Image type returned by computeImage.
Definition: Psf.h:83
Transform LSST spatial data, such as lsst::geom::Point2D and lsst::geom::SpherePoint,...
Definition: Transform.h:68
Describe the colour of a source.
Definition: Color.h:26
_view_t::x_iterator x_iterator
An iterator for traversing the pixels in a row.
Definition: ImageBase.h:133
Parameters to control convolution.
Definition: warpExposure.h:250
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
int getHeight() const noexcept
Definition: Box.h:188
Matrix const & getMatrix() const noexcept
An intermediate base class for Psfs that use an image representation.
Definition: ImagePsf.h:40
A Psf class that maps an arbitrary Psf through a coordinate transformation.
Definition: WarpedPsf.h:50
boost::shared_ptr< afw::detection::Psf const > _undistortedPsf
Definition: WarpedPsf.h:83
boost::shared_ptr< afw::geom::TransformPoint2ToPoint2 const > _distortion
Definition: WarpedPsf.h:84
virtual geom::Box2I doComputeBBox(geom::Point2D const &position, afw::image::Color const &color) const
Definition: WarpedPsf.cc:217
virtual boost::shared_ptr< afw::detection::Psf > clone() const
Polymorphic deep copy. Usually unnecessary, as Psfs are immutable.
Definition: WarpedPsf.cc:176
WarpedPsf(boost::shared_ptr< afw::detection::Psf const > undistortedPsf, boost::shared_ptr< afw::geom::TransformPoint2ToPoint2 const > distortion, boost::shared_ptr< afw::math::WarpingControl const > control)
Construct WarpedPsf from unwarped psf and distortion.
Definition: WarpedPsf.cc:138
virtual boost::shared_ptr< afw::detection::Psf > resized(int width, int height) const
Return a clone with specified kernel dimensions.
Definition: WarpedPsf.cc:180
virtual geom::Point2D getAveragePosition() const
Return the average of the positions of the stars that went into this Psf.
Definition: WarpedPsf.cc:172
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.
lsst::geom::AffineTransform linearizeTransform(TransformPoint2ToPoint2 const &original, lsst::geom::Point2D const &inPoint)
Approximate a Transform by its local linearization.
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
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 > Point2D
Definition: Point.h:324
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.