LSSTApplications  17.0+11,17.0+34,17.0+56,17.0+57,17.0+59,17.0+7,17.0-1-g377950a+33,17.0.1-1-g114240f+2,17.0.1-1-g4d4fbc4+28,17.0.1-1-g55520dc+49,17.0.1-1-g5f4ed7e+52,17.0.1-1-g6dd7d69+17,17.0.1-1-g8de6c91+11,17.0.1-1-gb9095d2+7,17.0.1-1-ge9fec5e+5,17.0.1-1-gf4e0155+55,17.0.1-1-gfc65f5f+50,17.0.1-1-gfc6fb1f+20,17.0.1-10-g87f9f3f+1,17.0.1-11-ge9de802+16,17.0.1-16-ga14f7d5c+4,17.0.1-17-gc79d625+1,17.0.1-17-gdae4c4a+8,17.0.1-2-g26618f5+29,17.0.1-2-g54f2ebc+9,17.0.1-2-gf403422+1,17.0.1-20-g2ca2f74+6,17.0.1-23-gf3eadeb7+1,17.0.1-3-g7e86b59+39,17.0.1-3-gb5ca14a,17.0.1-3-gd08d533+40,17.0.1-30-g596af8797,17.0.1-4-g59d126d+4,17.0.1-4-gc69c472+5,17.0.1-6-g5afd9b9+4,17.0.1-7-g35889ee+1,17.0.1-7-gc7c8782+18,17.0.1-9-gc4bbfb2+3,w.2019.22
LSSTDataManagementBasePackage
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 
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 {
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 
217 geom::Box2I WarpedPsf::doComputeBBox(geom::Point2D const &position, afw::image::Color const &color) const {
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
Extent< int, N > ceil(Extent< double, N > const &input) noexcept
Return the component-wise ceil (round towards more positive).
Definition: Extent.cc:118
lsst::geom::AffineTransform linearizeTransform(TransformPoint2ToPoint2 const &original, lsst::geom::Point2D const &inPoint)
Approximate a Transform by its local linearization.
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
An affine coordinate transformation consisting of a linear transformation and an offset.
#define CONST_PTR(...)
A shared pointer to a const object.
Definition: base.h:47
table::Key< int > b
_view_t::x_iterator x_iterator
An iterator for traversing the pixels in a row.
Definition: ImageBase.h:134
int y
Definition: SpanSet.cc:49
table::Key< int > a
Parameters to control convolution.
Definition: warpExposure.h:250
Point< double, 2 > Point2D
Definition: Point.h:324
STL class.
T min(T... args)
A base class for image defects.
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.
virtual boost::shared_ptr< afw::detection::Psf > resized(int width, int height) const
Return a clone with specified kernel dimensions.
Definition: WarpedPsf.cc:180
Reports errors in the logical structure of the program.
Definition: Runtime.h:46
table::Box2IKey bbox
Definition: Detector.cc:169
T max(T... args)
Point< int, 2 > Point2I
Definition: Point.h:321
LinearTransform const & getLinear() const noexcept
virtual geom::Point2D getAveragePosition() const
Return the average of the positions of the stars that went into this Psf.
Definition: WarpedPsf.cc:172
virtual boost::shared_ptr< afw::detection::Psf > clone() const
Polymorphic deep copy. Usually unnecessary, as Psfs are immutable.
Definition: WarpedPsf.cc:176
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
Extent< int, 2 > Extent2I
Definition: Extent.h:397
Reports invalid arguments.
Definition: Runtime.h:66
#define PTR(...)
Definition: base.h:41
An intermediate base class for Psfs that use an image representation.
Definition: ImagePsf.h:40
Extent< int, N > floor(Extent< double, N > const &input) noexcept
Return the component-wise floor (round towards more negative).
Definition: Extent.cc:109
Describe the colour of a source.
Definition: Color.h:26
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects...
A polymorphic base class for representing an image&#39;s Point Spread Function.
Definition: Psf.h:76
boost::shared_ptr< afw::detection::Psf const > _undistortedPsf
Definition: WarpedPsf.h:83
An integer coordinate rectangle.
Definition: Box.h:54
AffineTransform const inverted() const
Return the inverse transform.
boost::shared_ptr< afw::geom::TransformPoint2ToPoint2 const > _distortion
Definition: WarpedPsf.h:84
virtual boost::shared_ptr< afw::detection::Psf::Image > doComputeKernelImage(geom::Point2D const &position, afw::image::Color const &color) const
These virtual member functions are private, not protected, because we only want derived classes to im...
Definition: WarpedPsf.cc:188
image::Image< Pixel > Image
Image type returned by computeImage.
Definition: Psf.h:85
A Psf class that maps an arbitrary Psf through a coordinate transformation.
Definition: WarpedPsf.h:50
Transform LSST spatial data, such as lsst::geom::Point2D and lsst::geom::SpherePoint, using an AST mapping.
Definition: Transform.h:67
std::shared_ptr< TransformPoint2ToPoint2 > makeTransform(lsst::geom::AffineTransform const &affine)
Wrap an lsst::geom::AffineTransform as a Transform.