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
ChebyshevBoundedField.cc
Go to the documentation of this file.
1 // -*- LSST-C++ -*-
2 /*
3  * LSST Data Management System
4  * Copyright 2008-2014 LSST Corporation.
5  *
6  * This product includes software developed by the
7  * LSST Project (http://www.lsst.org/).
8  *
9  * This program is free software: you can redistribute it and/or modify
10  * it under the terms of the GNU General Public License as published by
11  * the Free Software Foundation, either version 3 of the License, or
12  * (at your option) any later version.
13  *
14  * This program is distributed in the hope that it will be useful,
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17  * GNU General Public License for more details.
18  *
19  * You should have received a copy of the LSST License Statement and
20  * the GNU General Public License along with this program. If not,
21  * see <http://www.lsstcorp.org/LegalNotices/>.
22  */
23 
24 #include "boost/make_shared.hpp"
25 
26 #include "ndarray/eigen.h"
34 
35 namespace lsst { namespace afw { namespace math {
36 
38  return detail::TrapezoidalPacker(*this).size;
39 }
40 
41 // ------------------ Constructors and helpers ---------------------------------------------------------------
42 
43 namespace {
44 
45 // Compute an affine transform that maps an arbitrary box to [-1,1]x[-1,1]
46 geom::AffineTransform makeChebyshevRangeTransform(afw::geom::Box2D const bbox) {
47  return geom::AffineTransform(
48  geom::LinearTransform::makeScaling(2.0 / bbox.getWidth(), 2.0 / bbox.getHeight()),
50  -(2.0*bbox.getCenterX()) / bbox.getWidth(),
51  -(2.0*bbox.getCenterY()) / bbox.getHeight()
52  )
53  );
54 }
55 
56 } // anonyomous
57 
59  afw::geom::Box2I const & bbox,
61 ) : BoundedField(bbox),
62  _toChebyshevRange(makeChebyshevRangeTransform(geom::Box2D(bbox))),
63  _coefficients(coefficients)
64 {}
65 
67  afw::geom::Box2I const & bbox
68 ) : BoundedField(bbox),
69  _toChebyshevRange(makeChebyshevRangeTransform(geom::Box2D(bbox)))
70 {}
71 
72 // ------------------ fit() and helpers ---------------------------------------------------------------------
73 
74 namespace {
75 
76 typedef ChebyshevBoundedField::Control Control;
77 typedef detail::TrapezoidalPacker Packer;
78 
79 // fill an array with 1-d Chebyshev functions of the 1st kind T(x), evaluated at the given point x
80 void evaluateBasis1d(
82  double x
83 ) {
84  int const n = t.getSize<0>();
85  if (n > 0) {
86  t[0] = 1.0;
87  }
88  if (n > 1) {
89  t[1] = x;
90  }
91  for (int i = 2; i < n; ++i) {
92  t[i] = 2.0*x*t[i-1] - t[i-2];
93  }
94 }
95 
96 // Create a matrix of 2-d Chebyshev functions evaluated at a set of positions, with
97 // Chebyshev order along columns and evaluation positions along rows. We pack the
98 // 2-d functions using the TrapezoidalPacker class, because we don't want any columns
99 // that correspond to coefficients that should be set to zero.
100 ndarray::Array<double,2,2> makeMatrix(
103  geom::AffineTransform const & toChebyshevRange,
104  Packer const & packer,
105  Control const & ctrl
106 ) {
107  int const nPoints = x.getSize<0>();
110  ndarray::Array<double,2,2> out = ndarray::allocate(nPoints, packer.size);
111  // Loop over x and y together, computing T_i(x) and T_j(y) arrays for each point,
112  // then packing them together.
113  for (int p = 0; p < nPoints; ++p) {
114  geom::Point2D sxy = toChebyshevRange(geom::Point2D(x[p], y[p]));
115  evaluateBasis1d(tx, sxy.getX());
116  evaluateBasis1d(ty, sxy.getY());
117  packer.pack(out[p], tx, ty); // this sets a row of out to the packed outer product of tx and ty
118  }
119  return out;
120 }
121 
122 // Create a matrix of 2-d Chebyshev functions evaluated on a grid of positions, with
123 // Chebyshev order along columns and evaluation positions along rows. We pack the
124 // 2-d functions using the TrapezoidalPacker class, because we don't want any columns
125 // that correspond to coefficients that should be set to zero.
126 ndarray::Array<double,2,2> makeMatrix(
127  geom::Box2I const & bbox,
128  geom::AffineTransform const & toChebyshevRange,
129  Packer const & packer,
130  Control const & ctrl
131 ) {
132  // Create a 2-d array that contains T_j(x) for each x value, with x values in rows and j in columns
133  ndarray::Array<double,2,2> tx = ndarray::allocate(bbox.getWidth(), packer.nx);
134  for (int x = bbox.getBeginX(), p = 0; p < bbox.getWidth(); ++p, ++x) {
135  evaluateBasis1d(
136  tx[p],
137  toChebyshevRange[geom::AffineTransform::XX]*x + toChebyshevRange[geom::AffineTransform::X]
138  );
139  }
140 
141  // Loop over y values, and at each point, compute T_i(y), then loop over x and multiply by the T_j(x)
142  // we already computed and stored above.
143  ndarray::Array<double,2,2> out = ndarray::allocate(bbox.getArea(),packer.size);
145  ndarray::Array<double,1,1> ty = ndarray::allocate(ctrl.orderY + 1);
146  for (int y = bbox.getBeginY(), i = 0; i < bbox.getHeight(); ++i, ++y) {
147  evaluateBasis1d(
148  ty,
149  toChebyshevRange[geom::AffineTransform::YY]*y + toChebyshevRange[geom::AffineTransform::Y]
150  );
151  for (int j = 0; j < bbox.getWidth(); ++j, ++outIter) {
152  // this sets a row of out to the packed outer product of tx and ty
153  packer.pack(*outIter, tx[j], ty);
154  }
155  }
156  return out;
157 }
158 
159 } // anonymous
160 
162  afw::geom::Box2I const & bbox,
163  ndarray::Array<double const,1> const & x,
164  ndarray::Array<double const,1> const & y,
165  ndarray::Array<double const,1> const & z,
166  Control const & ctrl
167 ) {
168  // Initialize the result object, so we can make use of the AffineTransform it builds
170  // This packer object knows how to map the 2-d Chebyshev functions onto a 1-d array,
171  // using only those that the control says should have nonzero coefficients.
172  Packer const packer(ctrl);
173  // Create a "design matrix" for the linear least squares problem (A in min||Ax-b||)
174  ndarray::Array<double,2,2> matrix = makeMatrix(x, y, result->_toChebyshevRange, packer, ctrl);
175  // Solve the linear least squares problem.
177  // Unpack the solution into a 2-d matrix, with zeros for values we didn't fit.
178  result->_coefficients = packer.unpack(lstsq.getSolution());
179  return result;
180 }
181 
183  afw::geom::Box2I const & bbox,
184  ndarray::Array<double const,1> const & x,
185  ndarray::Array<double const,1> const & y,
186  ndarray::Array<double const,1> const & z,
187  ndarray::Array<double const,1> const & w,
188  Control const & ctrl
189 ) {
190  // Initialize the result object, so we can make use of the AffineTransform it builds
192  // This packer object knows how to map the 2-d Chebyshev functions onto a 1-d array,
193  // using only those that the control says should have nonzero coefficients.
194  Packer const packer(ctrl);
195  // Create a "design matrix" for the linear least squares problem ('A' in min||Ax-b||)
196  ndarray::Array<double,2,2> matrix = makeMatrix(x, y, result->_toChebyshevRange, packer, ctrl);
197  // We want to do weighted least squares, so we multiply both the data vector 'b' and the
198  // matrix 'A' by the weights.
199  matrix.asEigen<Eigen::ArrayXpr>().colwise() *= w.asEigen<Eigen::ArrayXpr>();
201  wz.asEigen<Eigen::ArrayXpr>() *= w.asEigen<Eigen::ArrayXpr>();
202  // Solve the linear least squares problem.
204  // Unpack the solution into a 2-d matrix, with zeros for values we didn't fit.
205  result->_coefficients = packer.unpack(lstsq.getSolution());
206  return result;
207 }
208 
209 template <typename T>
211  image::Image<T> const & img,
212  Control const & ctrl
213 ) {
214  // Initialize the result object, so we can make use of the AffineTransform it builds
215  geom::Box2I bbox = img.getBBox(image::PARENT);
217  // This packer object knows how to map the 2-d Chebyshev functions onto a 1-d array,
218  // using only those that the control says should have nonzero coefficients.
219  Packer const packer(ctrl);
220  ndarray::Array<double,2,2> matrix = makeMatrix(bbox, result->_toChebyshevRange, packer, ctrl);
221  // Flatten the data image into a 1-d vector.
222  ndarray::Array<double,2,2> imgCopy = ndarray::allocate(img.getArray().getShape());
223  imgCopy.deep() = img.getArray();
224  ndarray::Array<double const,1,1> z = ndarray::flatten<1>(imgCopy);
225  // Solve the linear least squares problem.
227  // Unpack the solution into a 2-d matrix, with zeros for values we didn't fit.
228  result->_coefficients = packer.unpack(lstsq.getSolution());
229  return result;
230 }
231 
232 // ------------------ modifier factories ---------------------------------------------------------------
233 
235  if (ctrl.orderX >= _coefficients.getSize<1>()) {
236  throw LSST_EXCEPT(
237  pex::exceptions::LengthError,
238  (boost::format("New x order (%d) exceeds old x order (%d)")
239  % ctrl.orderX % (_coefficients.getSize<1>() - 1)).str()
240  );
241  }
242  if (ctrl.orderY >= _coefficients.getSize<0>()) {
243  throw LSST_EXCEPT(
244  pex::exceptions::LengthError,
245  (boost::format("New y order (%d) exceeds old y order (%d)")
246  % ctrl.orderY % (_coefficients.getSize<0>() - 1)).str()
247  );
248  }
249  ndarray::Array<double,2,2> coefficients = ndarray::allocate(ctrl.orderY + 1, ctrl.orderX + 1);
250  coefficients.deep()
251  = _coefficients[ndarray::view(0, ctrl.orderY + 1)(0, ctrl.orderX + 1)];
252  if (ctrl.triangular) {
253  Packer packer(ctrl);
254  ndarray::Array<double,1,1> packed = ndarray::allocate(packer.size);
255  packer.pack(packed, coefficients);
256  packer.unpack(coefficients, packed);
257  }
258  return boost::make_shared<ChebyshevBoundedField>(getBBox(), coefficients);
259 }
260 
261 PTR(ChebyshevBoundedField) ChebyshevBoundedField::relocate(geom::Box2I const & bbox) const {
262  return boost::make_shared<ChebyshevBoundedField>(bbox, _coefficients);
263 }
264 
265 
266 // ------------------ evaluate() and helpers ---------------------------------------------------------------
267 
268 namespace {
269 
270 // To evaluate a 1-d Chebyshev function without needing to have workspace, we use the
271 // Clenshaw algorith, which is like going through the recurrence relation in reverse.
272 // The CoeffGetter argument g is something that behaves like an array, providing access
273 // to the coefficients.
274 template <typename CoeffGetter>
275 double evaluateFunction1d(
276  CoeffGetter g,
277  double x,
278  int size
279 ) {
280  double b_kp2=0.0, b_kp1=0.0;
281  for (int k = (size - 1); k > 0; --k) {
282  double b_k = g[k] + 2*x*b_kp1 - b_kp2;
283  b_kp2 = b_kp1;
284  b_kp1 = b_k;
285  }
286  return g[0] + x*b_kp1 - b_kp2;
287 }
288 
289 // This class imitates a 1-d array, by running evaluateFunction1d on a nested dimension;
290 // this lets us reuse the logic in evaluateFunction1d for both dimensions. Essentially,
291 // we run evaluateFunction1d on a column of coefficients to evaluate T_j(x), then pass
292 // the result of that to evaluateFunction1d with the results as the "coefficients" associated
293 // with the T_j(y) functions.
294 struct RecursionArrayImitator {
295 
296  double operator[](int i) const {
297  return evaluateFunction1d(coefficients[i], x, coefficients.getSize<1>());
298  }
299 
300  RecursionArrayImitator(ndarray::Array<double const,2,2> const & coefficients_, double x_) :
301  coefficients(coefficients_), x(x_)
302  {}
303 
305  double x;
306 };
307 
308 } // anonymous
309 
310 double ChebyshevBoundedField::evaluate(geom::Point2D const & position) const {
311  geom::Point2D p = _toChebyshevRange(position);
312  return evaluateFunction1d(RecursionArrayImitator(_coefficients, p.getX()),
313  p.getY(), _coefficients.getSize<0>());
314 }
315 
316 // ------------------ persistence ---------------------------------------------------------------------------
317 
318 namespace {
319 
320 struct PersistenceHelper {
321  table::Schema schema;
322  table::Key<int> orderX;
323  table::PointKey<int> bboxMin;
324  table::PointKey<int> bboxMax;
325  table::Key< table::Array<double> > coefficients;
326 
327  PersistenceHelper(int nx, int ny) :
328  schema(),
329  orderX(schema.addField<int>("order_x", "maximum Chebyshev function order in x")),
330  bboxMin(table::PointKey<int>::addFields(
331  schema, "bbox_min", "lower-left corner of bounding box", "pixels")),
332  bboxMax(table::PointKey<int>::addFields(
333  schema, "bbox_max", "upper-right corner of bounding box", "pixels")),
334  coefficients(
335  schema.addField< table::Array<double> >(
336  "coefficients", "Chebyshev function coefficients, ordered by y then x", nx*ny
337  )
338  )
339  {}
340 
341  PersistenceHelper(table::Schema const & s) :
342  schema(s),
343  orderX(s["order_x"]),
344  bboxMin(s["bbox_min"]),
345  bboxMax(s["bbox_max"]),
346  coefficients(s["coefficients"])
347  {}
348 
349 };
350 
351 class ChebyshevBoundedFieldFactory : public table::io::PersistableFactory {
352 public:
353 
354  virtual PTR(table::io::Persistable)
355  read(InputArchive const & archive, CatalogVector const & catalogs) const {
356  LSST_ARCHIVE_ASSERT(catalogs.size() == 1u);
357  LSST_ARCHIVE_ASSERT(catalogs.front().size() == 1u);
358  table::BaseRecord const & record = catalogs.front().front();
359  PersistenceHelper const keys(record.getSchema());
360  geom::Box2I bbox(record.get(keys.bboxMin), record.get(keys.bboxMax));
361  int nx = record.get(keys.orderX) + 1;
362  int ny = keys.coefficients.getSize() / nx;
363  LSST_ARCHIVE_ASSERT(nx * ny == keys.coefficients.getSize());
365  ndarray::flatten<1>(coefficients) = record.get(keys.coefficients);
366  return boost::make_shared<ChebyshevBoundedField>(bbox, coefficients);
367  }
368 
369  ChebyshevBoundedFieldFactory(std::string const & name) : afw::table::io::PersistableFactory(name) {}
370 
371 };
372 
373 std::string getChebyshevBoundedFieldPersistenceName() { return "ChebyshevBoundedField"; }
374 
375 ChebyshevBoundedFieldFactory registration(getChebyshevBoundedFieldPersistenceName());
376 
377 } // anonymous
378 
380  return getChebyshevBoundedFieldPersistenceName();
381 }
382 
384  return "lsst.afw.math";
385 }
386 
388  PersistenceHelper const keys(_coefficients.getSize<1>(), _coefficients.getSize<0>());
389  table::BaseCatalog catalog = handle.makeCatalog(keys.schema);
390  PTR(table::BaseRecord) record = catalog.addNew();
391  record->set(keys.orderX, _coefficients.getSize<1>() - 1);
392  record->set(keys.bboxMin, getBBox().getMin());
393  record->set(keys.bboxMax, getBBox().getMax());
394  (*record)[keys.coefficients].deep() = ndarray::flatten<1>(_coefficients);
395  handle.saveCatalog(catalog);
396 }
397 
398 PTR(BoundedField) ChebyshevBoundedField::operator*(double const scale) const {
399  return boost::make_shared<ChebyshevBoundedField>(getBBox(), ndarray::copy(getCoefficients()*scale));
400 }
401 
402 // ------------------ explicit instantiation ----------------------------------------------------------------
403 
404 #ifndef DOXYGEN
405 
406 #define INSTANTIATE(T) \
407  template PTR(ChebyshevBoundedField) ChebyshevBoundedField::fit( \
408  image::Image<T> const & image, \
409  Control const & ctrl \
410  )
411 
412 INSTANTIATE(float);
413 INSTANTIATE(double);
414 
415 #endif
416 
417 }}} // namespace lsst::afw::math
int y
Extent< int, N > truncate(Extent< double, N > const &input)
View< boost::fusion::vector1< index::Full > > view()
Start a view definition that includes the entire first dimension.
Definition: views.h:131
int computeSize() const
Return the number of nonzero coefficients in the Chebyshev function defined by this object...
static LeastSquares fromDesignMatrix(ndarray::Array< T1, 2, C1 > const &design, ndarray::Array< T2, 1, C2 > const &data, Factorization factorization=NORMAL_EIGENSYSTEM)
Initialize from the design matrix and data vector given as ndarrays.
Definition: LeastSquares.h:99
#define PTR(...)
Definition: base.h:41
table::Key< std::string > name
Definition: ApCorrMap.cc:71
BaseCatalog makeCatalog(Schema const &schema)
Return a new, empty catalog with the given schema.
Eigen matrix objects that present a view into an ndarray::Array.
An object passed to Persistable::write to allow it to persist itself.
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
Definition: Persistable.h:47
virtual std::string getPersistenceName() const
Return the unique name used to persist this object and look up its factory.
int getBeginX() const
Definition: Box.h:139
ChebyshevBoundedField(afw::geom::Box2I const &bbox, ndarray::Array< double const, 2, 2 > const &coefficients)
Initialize the field from its bounding box an coefficients.
A custom container class for records, based on std::vector.
Definition: Catalog.h:94
A control object used when fitting ChebyshevBoundedField to data (see ChebyshevBoundedField::fit) ...
SelectEigenView< T >::Type copy(Eigen::EigenBase< T > const &other)
Copy an arbitrary Eigen expression into a new EigenView.
Definition: eigen.h:390
EigenView< Element, ND::value, RMC::value, XprKind, Rows, Cols > asEigen() const
double getHeight() const
Definition: Box.h:360
Solver for linear least-squares problems.
Definition: LeastSquares.h:65
int getSize() const
Return the size of a specific dimension.
Definition: ArrayBase.h:126
boost::shared_ptr< RecordT > addNew()
Create a new record, add it to the end of the catalog, and return a pointer to it.
Definition: Catalog.h:469
An integer coordinate rectangle.
Definition: Box.h:53
table::Key< table::Array< Kernel::Pixel > > image
Definition: FixedKernel.cc:117
double getCenterX() const
Definition: Box.h:374
int getHeight() const
Definition: Box.h:155
double w
Definition: CoaddPsf.cc:57
An affine coordinate transformation consisting of a linear transformation and an offset.
virtual void write(OutputArchiveHandle &handle) const
Write the object to one or more catalogs.
double getWidth() const
Definition: Box.h:359
static LinearTransform makeScaling(double s)
geom::Box2I getBBox() const
Definition: BoundedField.h:95
virtual double evaluate(geom::Point2D const &position) const
int getArea() const
Definition: Box.h:156
tbl::Schema schema
int x
detail::SimpleInitializer< N > allocate(Vector< int, N > const &shape)
Create an expression that allocates uninitialized memory for an array.
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
tbl::PointKey< int > bboxMax
#define INSTANTIATE(T)
Deep const deep() const
Return an ArrayRef view to this.
Definition: ArrayBase.h:178
Base class for all records.
Definition: BaseRecord.h:27
A BoundedField based on 2-d Chebyshev polynomials of the first kind.
Use the normal equations with a symmetric Eigensystem decomposition.
Definition: LeastSquares.h:71
An abstract base class for 2-d functions defined on an integer bounding boxes.
Definition: BoundedField.h:49
ndarray::Array< double const, 2, 2 > _coefficients
double getCenterY() const
Definition: Box.h:375
int getBeginY() const
Definition: Box.h:140
ExpressionTraits< Derived >::Iterator Iterator
Nested expression or element iterator.
A floating-point coordinate rectangle geometry.
Definition: Box.h:271
table::Key< int > orderX
ndarray::Array< double const, 2, 2 > coefficients
virtual std::string getPythonModule() const
Return the fully-qualified Python module that should be imported to guarantee that its factory is reg...
void saveCatalog(BaseCatalog const &catalog)
Save a catalog in the archive.
ndarray::Array< double const, 1, 1 > getSolution()
Return the vector solution to the least squares problem.
int getWidth() const
Definition: Box.h:154
Iterator begin() const
Return an Iterator to the beginning of the array.
Definition: ArrayBase.h:99
tbl::PointKey< int > bboxMin