LSST Applications g2079a07aa2+86d27d4dc4,g2305ad1205+a659bff248,g2bbee38e9b+3c60f8fe34,g337abbeb29+3c60f8fe34,g33d1c0ed96+3c60f8fe34,g3502564af9+d77d6d1350,g3a166c0a6a+3c60f8fe34,g487adcacf7+25d9892218,g4be5004598+d77d6d1350,g50ff169b8f+96c6868917,g52b1c1532d+585e252eca,g591dd9f2cf+4d81263f9a,g5cd07815a0+980d2b1c3b,g607f77f49a+d77d6d1350,g858d7b2824+d77d6d1350,g88963caddf+83e433e629,g99cad8db69+a4d3c48eeb,g9ddcbc5298+9a081db1e4,ga1e77700b3+bcf1af89ad,ga57fefb910+9a39d7b2d7,gae0086650b+585e252eca,gb065fddaf9+4f9fd82a2c,gb0e22166c9+60f28cb32d,gb363559e06+d84b1d3d07,gb3b7280ab2+4563d032e1,gb4b16eec92+babe958938,gba4ed39666+c2a2e4ac27,gbb8dafda3b+ed6854b564,gc120e1dc64+b72d212f87,gc28159a63d+3c60f8fe34,gc3e9b769f7+921dbcd359,gcf0d15dbbd+9a39d7b2d7,gdaeeff99f8+f9a426f77a,gddc38dedce+585e252eca,ge79ae78c31+3c60f8fe34,w.2024.21
LSST Data Management Base Package
Loading...
Searching...
No Matches
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"
31
32namespace lsst {
33namespace afw {
34
37
38namespace math {
39
40ndarray::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
49double BoundedField::integrate() const { throw LSST_EXCEPT(pex::exceptions::LogicError, "Not Implemented"); }
50
51double BoundedField::mean() const { throw LSST_EXCEPT(pex::exceptions::LogicError, "Not Implemented"); }
52
53namespace {
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
61struct Assign {
62 template <typename T, typename U>
63 void operator()(T &out, U a) const {
64 out = a;
65 }
66};
67
68struct 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
79struct Multiply {
80 template <typename T, typename U>
81 void operator()(T &out, U a) const {
82 out *= a;
83 }
84};
85
86struct 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.
98class Interpolator {
99public:
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),
126 _z00(std::numeric_limits<double>::quiet_NaN()),
127 _z01(std::numeric_limits<double>::quiet_NaN()),
128 _z10(std::numeric_limits<double>::quiet_NaN()),
129 _z11(std::numeric_limits<double>::quiet_NaN()) {}
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
151private:
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
221template <typename T, typename F>
222void 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
258template <typename T>
259void BoundedField::fillImage(image::Image<T> &img, bool overlapOnly, int xStep, int yStep) const {
260 applyToImage(*this, img, Assign(), overlapOnly, xStep, yStep);
261}
262
263template <typename T>
264void 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
269template <typename T>
270void BoundedField::multiplyImage(image::Image<T> &img, bool overlapOnly, int xStep, int yStep) const {
271 applyToImage(*this, img, Multiply(), overlapOnly, xStep, yStep);
272}
273
274template <typename T>
275void 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
285INSTANTIATE(float);
286INSTANTIATE(double);
287
288} // namespace math
289} // namespace afw
290} // namespace lsst
int const step
int min
double scaleBy
int end
int max
#define INSTANTIATE(FROMSYS, TOSYS)
Definition Detector.cc:509
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
int y
Definition SpanSet.cc:48
table::Key< int > a
int getX0() const
Return the image's column-origin.
Definition ImageBase.h:306
int getWidth() const
Return the number of columns in the image.
Definition ImageBase.h:294
lsst::geom::Box2I getBBox(ImageOrigin origin=PARENT) const
Definition ImageBase.h:445
int getY0() const
Return the image's row-origin.
Definition ImageBase.h:314
x_iterator x_at(int x, int y) const
Return an x_iterator to the point (x, y) in the image.
Definition ImageBase.h:407
A class to represent a 2-dimensional array of pixels.
Definition Image.h:51
Image subset(lsst::geom::Box2I const &bbox, ImageOrigin origin=PARENT) const
Return a subimage corresponding to the given box.
Definition Image.h:219
static std::shared_ptr< T > dynamicCast(std::shared_ptr< Persistable > const &ptr)
Dynamically cast a shared_ptr.
An integer coordinate rectangle.
Definition Box.h:55
int getBeginX() const noexcept
Definition Box.h:172
int getEndY() const noexcept
Definition Box.h:177
int getBeginY() const noexcept
Definition Box.h:173
int getMaxX() const noexcept
Definition Box.h:161
int getMaxY() const noexcept
Definition Box.h:162
int getEndX() const noexcept
Definition Box.h:176
Reports errors in the logical structure of the program.
Definition Runtime.h:46
Eigen::Array< int, 4, 1 > Bounds
T iota(T... args)
STL namespace.