LSSTApplications  18.1.0
LSSTDataManagementBasePackage
BoundedField.cc
Go to the documentation of this file.
1 // -*- LSST-C++ -*-
2 
3 /*
4  * LSST Data Management System
5  * Copyright 2008-2014 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 <numeric>
26 
27 #include "lsst/pex/exceptions.h"
29 #include "lsst/afw/table/io/Persistable.cc"
31 
32 namespace lsst {
33 namespace afw {
34 
37 
38 namespace math {
39 
40 ndarray::Array<double, 1, 1> BoundedField::evaluate(ndarray::Array<double const, 1> const &x,
41  ndarray::Array<double const, 1> const &y) const {
42  ndarray::Array<double, 1, 1> out = ndarray::allocate(x.getSize<0>());
43  for (int i = 0, n = x.getSize<0>(); i < n; ++i) {
44  out[i] = evaluate(x[i], y[i]);
45  }
46  return out;
47 }
48 
49 double BoundedField::integrate() const { throw LSST_EXCEPT(pex::exceptions::LogicError, "Not Implemented"); }
50 
51 double BoundedField::mean() const { throw LSST_EXCEPT(pex::exceptions::LogicError, "Not Implemented"); }
52 
53 namespace {
54 
55 // We use these operator-based functors to implement the various image-modifying routines
56 // in BoundedField. I don't think this is necessarily the best way to add interoperability
57 // with images, but it seems like a reasonable point on the simplicity vs. featurefulness
58 // scale, and I don't have time to explore more complex solutions (expression templates!).
59 // Of course, in C++11, we could replace all these with lambdas.
60 
61 struct Assign {
62  template <typename T, typename U>
63  void operator()(T &out, U a) const {
64  out = a;
65  }
66 };
67 
68 struct ScaledAdd {
69  explicit ScaledAdd(double s) : scaleBy(s) {}
70 
71  template <typename T, typename U>
72  void operator()(T &out, U a) const {
73  out += scaleBy * a;
74  }
75 
76  double scaleBy;
77 };
78 
79 struct Multiply {
80  template <typename T, typename U>
81  void operator()(T &out, U a) const {
82  out *= a;
83  }
84 };
85 
86 struct Divide {
87  template <typename T, typename U>
88  void operator()(T &out, U a) const {
89  out /= a;
90  }
91 };
92 
93 // Helper class to do bilinear interpolation with no dynamic memory
94 // allocation. This means we evaluate the function multiple times at some
95 // points (when we move to a new interval in y, to be specific). Could
96 // probably be optimized for the cases where either step is 1, but that seems
97 // premature.
98 class Interpolator {
99 public:
100  // Description of a cell to interpolate in one dimension.
101  struct Bounds {
102  int const step; // step between min and max, except in the special last row/column
103  int min; // lower-bound of cell (coordinate of known value and one before first point to fill in)
104  int max; // upper-bound of cell (coordinate of known value)
105  int end; // upper-bound of cell (one after last point to fill in)
106 
107  // Construct from step only.
108  //
109  // Other variables are initialized (and re-initialized) by calls to reset().
110  explicit Bounds(int step_) : step(step_), min(0), max(0), end(0) {}
111 
112  // Reset all points (aside from the step) to the first cell in this dimension.
113  void reset(int min_) {
114  min = min_;
115  max = min_ + step;
116  end = min_ + step;
117  }
118  };
119 
120  // Construct an object to interpolate the given BoundedField on an evenly-spaced grid within a region.
121  Interpolator(BoundedField const *field, lsst::geom::Box2I const *region, int xStep, int yStep)
122  : _field(field),
123  _region(region),
124  _x(xStep),
125  _y(yStep),
130 
131  // Actually do the interpolation.
132  //
133  // This method iterates over rows of cells in y, calling _runRow()
134  // to iterate over cells in x.
135  template <typename T, typename F>
136  void run(image::Image<T> &img, F functor) {
137  _y.reset(_region->getBeginY());
138  while (_y.end < _region->getEndY()) {
139  _runRow(img, functor);
140  _y.min = _y.max;
141  _y.max += _y.step;
142  _y.end = _y.max;
143  }
144  { // special-case last iteration in y
145  _y.max = _region->getMaxY();
146  _y.end = _region->getEndY();
147  _runRow(img, functor);
148  }
149  }
150 
151 private:
152  // Process a row of cells, calling _runCell() on each one.
153  template <typename T, typename F>
154  void _runRow(image::Image<T> &img, F functor) {
155  _x.reset(_region->getBeginX());
156  _z00 = _field->evaluate(_x.min, _y.min);
157  _z01 = _field->evaluate(_x.min, _y.max);
158  while (_x.max < _region->getEndX()) {
159  _z10 = _field->evaluate(_x.max, _y.min);
160  _z11 = _field->evaluate(_x.max, _y.max);
161  _runCell(img, functor);
162  _x.min = _x.max;
163  _x.max += _x.step;
164  _x.end = _x.max;
165  _z00 = _z10;
166  _z01 = _z11;
167  }
168  { // special-case last iteration in x
169  _x.max = _region->getMaxX();
170  _x.end = _region->getEndX();
171  _z10 = _field->evaluate(_x.max, _y.min);
172  _z11 = _field->evaluate(_x.max, _y.max);
173  _runCell(img, functor);
174  }
175  }
176 
177  // Interpolate all points in a cell, which is defined as a rectangle for which
178  // the BoundedField has been evaluated at all four corners. The main complication
179  // comes from the need to special-case the unusually-sized final cells and final
180  // rows/columns within cells in each dimension.
181  template <typename T, typename F>
182  void _runCell(image::Image<T> const &img, F functor) const {
183  int dy = _y.max - _y.min;
184  int dx = _x.max - _x.min;
185  // First iteration of each of the for loops below has been
186  // split into an explicit block. This avoids a few instructions
187  // when no interpolation is necessary, but more importantly it
188  // allows us to handle dx==0 and dy==0 with no divide-by-zero
189  // problems.
190  { // y=_y.min
191  auto rowIter = img.x_at(_x.min - img.getX0(), _y.min - img.getY0());
192  { // x=_x.min
193  functor(*rowIter, _z00);
194  ++rowIter;
195  }
196  for (int x = _x.min + 1; x < _x.end; ++x, ++rowIter) {
197  functor(*rowIter, ((_x.max - x) * _z00 + (x - _x.min) * _z10) / dx);
198  }
199  }
200  for (int y = _y.min + 1; y < _y.end; ++y) {
201  auto rowIter = img.x_at(_x.min - img.getX0(), y - img.getY0());
202  double z0 = ((_y.max - y) * _z00 + (y - _y.min) * _z01) / dy;
203  double z1 = ((_y.max - y) * _z10 + (y - _y.min) * _z11) / dy;
204  { // x=_x.min
205  functor(*rowIter, z0);
206  ++rowIter;
207  }
208  for (int x = _x.min + 1; x < _x.end; ++x, ++rowIter) {
209  functor(*rowIter, ((_x.max - x) * z0 + (x - _x.min) * z1) / dx);
210  }
211  }
212  }
213 
214  BoundedField const *_field;
215  lsst::geom::Box2I const *_region;
216  Bounds _x;
217  Bounds _y;
218  double _z00, _z01, _z10, _z11;
219 };
220 
221 template <typename T, typename F>
222 void applyToImage(BoundedField const &field, image::Image<T> &img, F functor, bool overlapOnly, int xStep,
223  int yStep) {
224  lsst::geom::Box2I region(field.getBBox());
225  if (overlapOnly) {
226  region.clip(img.getBBox(image::PARENT));
227  } else if (region != img.getBBox(image::PARENT)) {
228  throw LSST_EXCEPT(pex::exceptions::RuntimeError,
229  "Image bounding box does not match field bounding box");
230  }
231 
232  if (yStep > 1 || xStep > 1) {
233  Interpolator interpolator(&field, &region, xStep, yStep);
234  interpolator.run(img, functor);
235  } else {
236  // We iterate over rows as a significant optimization for AST-backed bounded fields
237  // (it's also slightly faster for other bounded fields, too).
238  auto subImage = img.subset(region);
239  auto size = region.getWidth();
240  ndarray::Array<double, 1> xx = ndarray::allocate(ndarray::makeVector(size));
241  ndarray::Array<double, 1> yy = ndarray::allocate(ndarray::makeVector(size));
242  // y gets incremented each outer loop, x is always xMin->xMax
243  std::iota(xx.begin(), xx.end(), region.getBeginX());
244  auto outRowIter = subImage.getArray().begin();
245  for (int y = region.getBeginY(); y < region.getEndY(); ++y, ++outRowIter) {
246  yy.deep() = y; // don't need indexToPosition, as we're already working in the right box (region).
247  functor(*outRowIter, field.evaluate(xx, yy));
248  }
249  }
250 }
251 
252 } // namespace
253 
255  return *bf * scale;
256 }
257 
258 template <typename T>
259 void BoundedField::fillImage(image::Image<T> &img, bool overlapOnly, int xStep, int yStep) const {
260  applyToImage(*this, img, Assign(), overlapOnly, xStep, yStep);
261 }
262 
263 template <typename T>
264 void BoundedField::addToImage(image::Image<T> &img, double scaleBy, bool overlapOnly, int xStep,
265  int yStep) const {
266  applyToImage(*this, img, ScaledAdd(scaleBy), overlapOnly, xStep, yStep);
267 }
268 
269 template <typename T>
270 void BoundedField::multiplyImage(image::Image<T> &img, bool overlapOnly, int xStep, int yStep) const {
271  applyToImage(*this, img, Multiply(), overlapOnly, xStep, yStep);
272 }
273 
274 template <typename T>
275 void BoundedField::divideImage(image::Image<T> &img, bool overlapOnly, int xStep, int yStep) const {
276  applyToImage(*this, img, Divide(), overlapOnly, xStep, yStep);
277 }
278 
279 #define INSTANTIATE(T) \
280  template void BoundedField::fillImage(image::Image<T> &, bool, int, int) const; \
281  template void BoundedField::addToImage(image::Image<T> &, double, bool, int, int) const; \
282  template void BoundedField::multiplyImage(image::Image<T> &, bool, int, int) const; \
283  template void BoundedField::divideImage(image::Image<T> &, bool, int, int) const
284 
285 INSTANTIATE(float);
286 INSTANTIATE(double);
287 
288 } // namespace math
289 } // namespace afw
290 } // namespace lsst
static std::shared_ptr< T > dynamicCast(std::shared_ptr< Persistable > const &ptr)
Dynamically cast a shared_ptr.
Definition: Persistable.cc:18
std::shared_ptr< Image< PixelT > > operator*(Image< PixelT > const &img, ImageSlice< PixelT > const &slc)
Overload operator*()
Definition: ImageSlice.cc:104
def scale(algorithm, min, max=None, frame=None)
Definition: ds9.py:109
int y
Definition: SpanSet.cc:49
#define INSTANTIATE(T)
table::Key< int > a
T iota(T... args)
int getX0() const
Return the image&#39;s column-origin.
Definition: ImageBase.h:324
int getEndY() const noexcept
Definition: Box.h:164
Image subset(lsst::geom::Box2I const &bbox, ImageOrigin origin=PARENT) const
Return a subimage corresponding to the given box.
Definition: Image.h:227
int min
int getBeginY() const noexcept
Definition: Box.h:160
A base class for image defects.
lsst::geom::Box2I getBBox(ImageOrigin origin=PARENT) const
Definition: ImageBase.h:463
int const step
int getMaxY() const noexcept
Definition: Box.h:149
int getBeginX() const noexcept
Definition: Box.h:159
int max
Reports errors in the logical structure of the program.
Definition: Runtime.h:46
solver_t * s
int getMaxX() const noexcept
Definition: Box.h:148
double x
double scaleBy
Definition: BoundedField.cc:76
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
int getY0() const
Return the image&#39;s row-origin.
Definition: ImageBase.h:332
int getEndX() const noexcept
Definition: Box.h:163
void addToImage(boost::shared_ptr< afw::image::Image< double > > image, std::vector< boost::shared_ptr< afw::image::Image< double > >> const &imgVector, std::vector< double > const &weightVector)
Definition: CoaddPsf.cc:192
x_iterator x_at(int x, int y) const
Return an x_iterator to the point (x, y) in the image.
Definition: ImageBase.h:425
An integer coordinate rectangle.
Definition: Box.h:54
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:59
int end
UnaryFunctionT::result_type integrate(UnaryFunctionT func, typename UnaryFunctionT::argument_type const a, typename UnaryFunctionT::argument_type const b, double eps=1.0e-6)
The 1D integrator.
Definition: Integrate.h:886