LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
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 
25 #include "lsst/pex/exceptions.h"
28 #include "lsst/afw/image/Image.h"
29 
30 namespace lsst { namespace meas { namespace algorithms {
31 
32 namespace {
33 
34 inline double min4(double a, double b, double c, double d) {
35  return std::min(std::min(a,b), std::min(c,d));
36 }
37 
38 inline double max4(double a, double b, double c, double d) {
39  return std::max(std::max(a,b), std::max(c,d));
40 }
41 
42 // TODO: make this routine externally callable and more generic using templates
43 // (also useful in e.g. math/offsetImage.cc)
44 PTR(afw::detection::Psf::Image) zeroPadImage(afw::detection::Psf::Image const &im, int xPad, int yPad) {
45  int nx = im.getWidth();
46  int ny = im.getHeight();
47 
48  PTR(afw::detection::Psf::Image) out = boost::make_shared<afw::detection::Psf::Image>(nx+2*xPad, ny+2*yPad);
49  out->setXY0(im.getX0()-xPad, im.getY0()-yPad);
50 
51  afw::geom::Box2I box(afw::geom::Point2I(xPad,yPad), afw::geom::Extent2I(nx,ny));
52  PTR(afw::detection::Psf::Image) subimage = boost::make_shared<afw::detection::Psf::Image>(
53  *out, box, lsst::afw::image::LOCAL);
54  *subimage <<= im;
55 
56  return out;
57 }
58 
69 PTR(afw::detection::Psf::Image) warpAffine(
70  afw::detection::Psf::Image const &im, afw::geom::AffineTransform const &t,
71  afw::math::WarpingControl const &wc
72 ) {
73  static const int dst_padding = 0;
74 
75  afw::math::SeparableKernel const& kernel = *wc.getWarpingKernel();
76  afw::geom::Point2I const& center = kernel.getCtr();
77  int const xPad = std::max(center.getX(), kernel.getWidth() - center.getX());
78  int const yPad = std::min(center.getY(), kernel.getHeight() - center.getY());
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(
88  pex::exceptions::RangeError,
89  "Unexpectedly large transform passed to WarpedPsf"
90  );
91  }
92 
93  // min/max coordinate values in input image
94  int in_xlo = im.getX0();
95  int in_xhi = im.getX0() + im.getWidth() - 1;
96  int in_ylo = im.getY0();
97  int in_yhi = im.getY0() + im.getHeight() - 1;
98 
99  // corners of output image
100  afw::geom::Point2D c00 = t(afw::geom::Point2D(in_xlo,in_ylo));
101  afw::geom::Point2D c01 = t(afw::geom::Point2D(in_xlo,in_yhi));
102  afw::geom::Point2D c10 = t(afw::geom::Point2D(in_xhi,in_ylo));
103  afw::geom::Point2D c11 = t(afw::geom::Point2D(in_xhi,in_yhi));
104 
105  //
106  // bounding box for output image
107  //
108  int out_xlo = floor(min4(c00.getX(),c01.getX(),c10.getX(),c11.getX())) - dst_padding;
109  int out_ylo = floor(min4(c00.getY(),c01.getY(),c10.getY(),c11.getY())) - dst_padding;
110  int out_xhi = ceil(max4(c00.getX(),c01.getX(),c10.getX(),c11.getX())) + dst_padding;
111  int out_yhi = ceil(max4(c00.getY(),c01.getY(),c10.getY(),c11.getY())) + dst_padding;
112 
113  // allocate output image
115  = boost::make_shared<afw::detection::Psf::Image>(out_xhi-out_xlo+1, out_yhi-out_ylo+1);
116  ret->setXY0(afw::geom::Point2I(out_xlo,out_ylo));
117 
118  // zero-pad input image
119  PTR(afw::detection::Psf::Image) im_padded = zeroPadImage(im, xPad, yPad);
120 
121  // warp it!
122  afw::math::warpImage(*ret, *im_padded, t, wc, 0.0);
123  return ret;
124 }
125 
126 } // anonymous
127 
129  PTR(afw::detection::Psf const) undistortedPsf,
130  PTR(afw::geom::XYTransform const) distortion,
131  CONST_PTR(afw::math::WarpingControl) control
132  ) :
133  ImagePsf(false),
134  _undistortedPsf(undistortedPsf),
135  _distortion(distortion),
136  _warpingControl(control)
137 {
138  _init();
139 }
140 
142  PTR(afw::detection::Psf const) undistortedPsf,
143  PTR(afw::geom::XYTransform const) distortion,
144  std::string const& kernelName,
145  unsigned int cache
146  ) :
147  ImagePsf(false),
148  _undistortedPsf(undistortedPsf),
149  _distortion(distortion),
150  _warpingControl(new afw::math::WarpingControl(kernelName, "", cache))
151 {
152  _init();
153 }
154 
156 {
157  if (!_undistortedPsf) {
158  throw LSST_EXCEPT(
159  pex::exceptions::LogicError,
160  "Undistorted Psf passed to WarpedPsf must not be None/NULL"
161  );
162  }
163  if (!_distortion) {
164  throw LSST_EXCEPT(
165  pex::exceptions::LogicError,
166  "XYTransform passed to WarpedPsf must not be None/NULL"
167  );
168  }
169  if (!_warpingControl) {
170  throw LSST_EXCEPT(
171  pex::exceptions::LogicError,
172  "WarpingControl passed to WarpedPsf must not be None/NULL"
173  );
174  }
175 }
176 
178  return _distortion->forwardTransform(_undistortedPsf->getAveragePosition());
179 }
180 
182  return boost::make_shared<WarpedPsf>(_undistortedPsf->clone(), _distortion->clone(), _warpingControl);
183 }
184 
185 PTR(afw::detection::Psf::Image) WarpedPsf::doComputeKernelImage(
186  afw::geom::Point2D const & position, afw::image::Color const & color
187 ) const {
188  afw::geom::AffineTransform t = _distortion->linearizeReverseTransform(position);
189  afw::geom::Point2D tp = t(position);
190 
191  PTR(Image) im = _undistortedPsf->computeKernelImage(tp, color);
192 
193  // Go to the warped coordinate system with 'p' at the origin
194  PTR(afw::detection::Psf::Psf::Image) ret
195  = warpAffine(*im, afw::geom::AffineTransform(t.invert().getLinear()), *_warpingControl);
196 
197  double normFactor = 1.0;
198  //
199  // Normalize the output image to sum 1
200  // FIXME defining a member function Image::getSum() would be convenient here and in other places
201  //
202  normFactor = 0.0;
203  for (int y = 0; y != ret->getHeight(); ++y) {
204  afw::detection::Psf::Image::x_iterator imEnd = ret->row_end(y);
205  for (afw::detection::Psf::Image::x_iterator imPtr = ret->row_begin(y); imPtr != imEnd; imPtr++) {
206  normFactor += *imPtr;
207  }
208  }
209  if (normFactor == 0.0) {
210  throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, "psf image has sum 0");
211  }
212  *ret /= normFactor;
213  return ret;
214 }
215 
216 }}} // namepace lsst::meas::algorithms
int y
virtual afw::geom::Point2D getAveragePosition() const
Return the average of the positions of the stars that went into this Psf.
Definition: WarpedPsf.cc:177
#define PTR(...)
Definition: base.h:41
AffineTransform const invert() const
boost::shared_ptr< afw::math::WarpingControl const > _warpingControl
Definition: WarpedPsf.h:92
#define CONST_PTR(...)
Definition: base.h:47
Point< int, 2 > Point2I
Definition: PSF.h:39
Extent< int, N > ceil(Extent< double, N > const &input)
int warpImage(DestImageT &destImage, lsst::afw::image::Wcs const &destWcs, SrcImageT const &srcImage, lsst::afw::image::Wcs 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. See also convenience function warpExposure() to warp an Ex...
table::Key< table::Array< Kernel::Pixel > > image
Definition: FixedKernel.cc:117
An affine coordinate transformation consisting of a linear transformation and an offset.
Support for 2-D images.
LinearTransform const & getLinear() const
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
Extent< int, N > floor(Extent< double, N > const &input)
Support for warping an image to a new WCS.
WarpedPsf(boost::shared_ptr< afw::detection::Psf const > undistortedPsf, boost::shared_ptr< afw::geom::XYTransform const > distortion, boost::shared_ptr< afw::math::WarpingControl const > control)
Construct WarpedPsf from unwarped psf and distortion.
Definition: WarpedPsf.cc:128
_view_t::x_iterator x_iterator
An iterator for traversing the pixels in a row.
Definition: Image.h:151
Virtual base class for 2D transforms.
Definition: XYTransform.h:48
afw::table::Key< double > b
Point< double, 2 > Point2D
Definition: Point.h:286
An intermediate base class for Psfs that use an image representation.
Definition: ImagePsf.h:37
A polymorphic base class for representing an image&#39;s Point Spread Function.
Definition: Psf.h:68
boost::shared_ptr< afw::detection::Psf const > _undistortedPsf
Definition: WarpedPsf.h:87
image::Image< Pixel > Image
Image type returned by computeImage.
Definition: Psf.h:79
Extent< int, 2 > Extent2I
Definition: Extent.h:355
A Psf class that maps an arbitrary Psf through a coordinate transformation.
Definition: WarpedPsf.h:49
boost::shared_ptr< afw::geom::XYTransform const > _distortion
Definition: WarpedPsf.h:88
Include files required for standard LSST Exception handling.