LSST Applications  21.0.0-147-g0e635eb1+1acddb5be5,22.0.0+052faf71bd,22.0.0+1ea9a8b2b2,22.0.0+6312710a6c,22.0.0+729191ecac,22.0.0+7589c3a021,22.0.0+9f079a9461,22.0.1-1-g7d6de66+b8044ec9de,22.0.1-1-g87000a6+536b1ee016,22.0.1-1-g8e32f31+6312710a6c,22.0.1-10-gd060f87+016f7cdc03,22.0.1-12-g9c3108e+df145f6f68,22.0.1-16-g314fa6d+c825727ab8,22.0.1-19-g93a5c75+d23f2fb6d8,22.0.1-19-gb93eaa13+aab3ef7709,22.0.1-2-g8ef0a89+b8044ec9de,22.0.1-2-g92698f7+9f079a9461,22.0.1-2-ga9b0f51+052faf71bd,22.0.1-2-gac51dbf+052faf71bd,22.0.1-2-gb66926d+6312710a6c,22.0.1-2-gcb770ba+09e3807989,22.0.1-20-g32debb5+b8044ec9de,22.0.1-23-gc2439a9a+fb0756638e,22.0.1-3-g496fd5d+09117f784f,22.0.1-3-g59f966b+1e6ba2c031,22.0.1-3-g849a1b8+f8b568069f,22.0.1-3-gaaec9c0+c5c846a8b1,22.0.1-32-g5ddfab5d3+60ce4897b0,22.0.1-4-g037fbe1+64e601228d,22.0.1-4-g8623105+b8044ec9de,22.0.1-5-g096abc9+d18c45d440,22.0.1-5-g15c806e+57f5c03693,22.0.1-7-gba73697+57f5c03693,master-g6e05de7fdc+c1283a92b8,master-g72cdda8301+729191ecac,w.2021.39
LSST Data Management Base Package
FunctionLibrary.h
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 #ifndef LSST_AFW_MATH_FUNCTIONLIBRARY_H
26 #define LSST_AFW_MATH_FUNCTIONLIBRARY_H
27 /*
28  * Define a collection of useful Functions.
29  */
30 #include <algorithm>
31 #include <cmath>
32 
33 #include "lsst/geom.h"
34 #include "lsst/afw/math/Function.h"
35 #include "lsst/geom/Angle.h"
36 
37 namespace lsst {
38 namespace afw {
39 namespace math {
40 
51 template <typename ReturnT>
52 class IntegerDeltaFunction1 : public Function1<ReturnT> {
53 public:
57  explicit IntegerDeltaFunction1(double xo) : Function1<ReturnT>(0), _xo(xo) {}
58 
63  ~IntegerDeltaFunction1() noexcept override = default;
64 
65  std::shared_ptr<Function1<ReturnT>> clone() const override {
67  }
68 
69  ReturnT operator()(double x) const noexcept(IS_NOTHROW_INIT<ReturnT>) override {
70  return static_cast<ReturnT>(x == _xo);
71  }
72 
73  std::string toString(std::string const& prefix = "") const override {
75  os << "IntegerDeltaFunction1 [" << _xo << "]: ";
76  os << Function1<ReturnT>::toString(prefix);
77  return os.str();
78  };
79 
80 private:
81  double _xo;
82 
83 protected:
84  /* Default constructor: intended only for serialization */
85  explicit IntegerDeltaFunction1() : Function1<ReturnT>(0), _xo(0.0) {}
86 };
87 
98 template <typename ReturnT>
99 class IntegerDeltaFunction2 : public Function2<ReturnT> {
100 public:
104  explicit IntegerDeltaFunction2(double xo, double yo) : Function2<ReturnT>(0), _xo(xo), _yo(yo) {}
105 
110  ~IntegerDeltaFunction2() noexcept override = default;
111 
112  std::shared_ptr<Function2<ReturnT>> clone() const override {
114  }
115 
116  ReturnT operator()(double x, double y) const noexcept(IS_NOTHROW_INIT<ReturnT>) override {
117  return static_cast<ReturnT>((x == _xo) && (y == _yo));
118  }
119 
120  std::string toString(std::string const& prefix) const override {
122  os << "IntegerDeltaFunction2 [" << _xo << ", " << _yo << "]: ";
123  os << Function2<ReturnT>::toString(prefix);
124  return os.str();
125  }
126 
127 private:
128  double _xo;
129  double _yo;
130 
131 protected:
132  /* Default constructor: intended only for serialization */
133  explicit IntegerDeltaFunction2() : Function2<ReturnT>(0), _xo(0.0), _yo(0.0) {}
134 };
135 
146 template <typename ReturnT>
147 class GaussianFunction1 : public Function1<ReturnT> {
148 public:
152  explicit GaussianFunction1(double sigma)
153  : Function1<ReturnT>(1), _multFac(1.0 / std::sqrt(lsst::geom::TWOPI)) {
154  this->_params[0] = sigma;
155  }
160  ~GaussianFunction1() noexcept override = default;
161 
162  std::shared_ptr<Function1<ReturnT>> clone() const override {
164  }
165 
166  ReturnT operator()(double x) const noexcept(IS_NOTHROW_INIT<ReturnT>) override {
167  return static_cast<ReturnT>((_multFac / this->_params[0]) *
168  std::exp(-(x * x) / (2.0 * this->_params[0] * this->_params[0])));
169  }
170 
171  std::string toString(std::string const& prefix) const override {
173  os << "GaussianFunction1 [" << _multFac << "]: ";
174  os << Function1<ReturnT>::toString(prefix);
175  return os.str();
176  }
177 
178 private:
179  const double _multFac;
180 
181 protected:
182  /* Default constructor: intended only for serialization */
183  explicit GaussianFunction1() : Function1<ReturnT>(1), _multFac(1.0 / std::sqrt(lsst::geom::TWOPI)) {}
184 };
185 
200 template <typename ReturnT>
201 class GaussianFunction2 : public Function2<ReturnT> {
202 public:
206  explicit GaussianFunction2(double sigma1,
207  double sigma2,
208  double angle = 0.0)
209  : Function2<ReturnT>(3), _multFac(1.0 / (lsst::geom::TWOPI)) {
210  this->_params[0] = sigma1;
211  this->_params[1] = sigma2;
212  this->_params[2] = angle;
213  _updateCache();
214  }
215 
220  ~GaussianFunction2() noexcept override = default;
221 
222  std::shared_ptr<Function2<ReturnT>> clone() const override {
224  new GaussianFunction2(this->_params[0], this->_params[1], this->_params[2]));
225  }
226 
227  ReturnT operator()(double x, double y) const noexcept(IS_NOTHROW_INIT<ReturnT>) override {
228  if (_angle != this->_params[2]) {
229  _updateCache();
230  }
231  double pos1 = (_cosAngle * x) + (_sinAngle * y);
232  double pos2 = (-_sinAngle * x) + (_cosAngle * y);
233  return static_cast<ReturnT>((_multFac / (this->_params[0] * this->_params[1])) *
234  std::exp(-((pos1 * pos1) / (2.0 * this->_params[0] * this->_params[0])) -
235  ((pos2 * pos2) / (2.0 * this->_params[1] * this->_params[1]))));
236  }
237 
238  std::string toString(std::string const& prefix) const override {
240  os << "GaussianFunction2: ";
241  os << Function2<ReturnT>::toString(prefix);
242  return os.str();
243  }
244 
245  bool isPersistable() const noexcept override { return true; }
246 
247 protected:
248  std::string getPersistenceName() const override;
249 
250  void write(afw::table::io::OutputArchiveHandle& handle) const override;
251 
252 private:
270  void _updateCache() const noexcept {
271  _angle = this->_params[2];
272  _sinAngle = std::sin(_angle);
273  _cosAngle = std::cos(_angle);
274  }
275  const double _multFac;
276  mutable double _angle;
277  mutable double _sinAngle;
278  mutable double _cosAngle;
279 
280 protected:
281  /* Default constructor: intended only for serialization */
282  explicit GaussianFunction2()
283  : Function2<ReturnT>(3),
284  _multFac(1.0 / (lsst::geom::TWOPI)),
285  _angle(0.0),
286  _sinAngle(0.0),
287  _cosAngle(1.0) {}
288 };
289 
304 template <typename ReturnT>
305 class DoubleGaussianFunction2 : public Function2<ReturnT> {
306 public:
311  double sigma1,
312  double sigma2 = 0,
313  double ampl2 = 0)
314  : Function2<ReturnT>(3), _multFac(1.0 / (lsst::geom::TWOPI)) {
315  this->_params[0] = sigma1;
316  this->_params[1] = sigma2;
317  this->_params[2] = ampl2;
318  }
319 
324  ~DoubleGaussianFunction2() noexcept override = default;
325 
326  std::shared_ptr<Function2<ReturnT>> clone() const override {
328  new DoubleGaussianFunction2(this->_params[0], this->_params[1], this->_params[2]));
329  }
330 
331  ReturnT operator()(double x, double y) const noexcept(IS_NOTHROW_INIT<ReturnT>) override {
332  double radSq = (x * x) + (y * y);
333  double sigma1Sq = this->_params[0] * this->_params[0];
334  double sigma2Sq = this->_params[1] * this->_params[1];
335  double b = this->_params[2];
336  return static_cast<ReturnT>(
337  (_multFac / (sigma1Sq + (b * sigma2Sq))) *
338  (std::exp(-radSq / (2.0 * sigma1Sq)) + (b * std::exp(-radSq / (2.0 * sigma2Sq)))));
339  }
340 
341  std::string toString(std::string const& prefix) const override {
343  os << "DoubleGaussianFunction2 [" << _multFac << "]: ";
344  os << Function2<ReturnT>::toString(prefix);
345  return os.str();
346  }
347 
348  bool isPersistable() const noexcept override { return true; }
349 
350 protected:
351  std::string getPersistenceName() const override;
352 
353  void write(afw::table::io::OutputArchiveHandle& handle) const override;
354 
355 private:
356  const double _multFac;
357 
358 protected:
359  /* Default constructor: intended only for serialization */
360  explicit DoubleGaussianFunction2() : Function2<ReturnT>(3), _multFac(1.0 / (lsst::geom::TWOPI)) {}
361 };
362 
370 template <typename ReturnT>
371 class PolynomialFunction1 : public Function1<ReturnT> {
372 public:
378  explicit PolynomialFunction1(unsigned int order)
379  : Function1<ReturnT>(order + 1) {}
380 
389  : Function1<ReturnT>(params) {
390  if (params.size() < 1) {
392  "PolynomialFunction1 called with empty vector");
393  }
394  }
395 
400  ~PolynomialFunction1() noexcept override = default;
401 
402  std::shared_ptr<Function1<ReturnT>> clone() const override {
404  }
405 
406  bool isLinearCombination() const noexcept override { return true; };
407 
408  ReturnT operator()(double x) const noexcept(IS_NOTHROW_INIT<ReturnT>) override {
409  int const order = static_cast<int>(this->_params.size()) - 1;
410  double retVal = this->_params[order];
411  for (int ii = order - 1; ii >= 0; --ii) {
412  retVal = (retVal * x) + this->_params[ii];
413  }
414  return static_cast<ReturnT>(retVal);
415  }
416 
420  unsigned int getOrder() const noexcept { return this->getNParameters() - 1; };
421 
422  std::string toString(std::string const& prefix) const override {
424  os << "PolynomialFunction1 []: ";
425  os << Function1<ReturnT>::toString(prefix);
426  return os.str();
427  }
428 
429 protected:
430  /* Default constructor: intended only for serialization */
431  explicit PolynomialFunction1() : Function1<ReturnT>(1) {}
432 };
433 
448 template <typename ReturnT>
450 public:
458  explicit PolynomialFunction2(unsigned int order)
459  : BasePolynomialFunction2<ReturnT>(order), _oldY(0), _xCoeffs(this->_order + 1) {}
460 
472  std::vector<double> params)
474  : BasePolynomialFunction2<ReturnT>(params), _oldY(0), _xCoeffs(this->_order + 1) {}
475 
480  ~PolynomialFunction2() noexcept override = default;
481 
482  std::shared_ptr<Function2<ReturnT>> clone() const override {
484  }
485 
486  ReturnT operator()(double x, double y) const noexcept(IS_NOTHROW_INIT<ReturnT>) override {
487  /* Solve as follows:
488  - f(x,y) = Cx0 + Cx1 x + Cx2 x^2 + Cx3 x^3 + ...
489  where:
490  Cx0 = P0 + P2 y + P5 y^2 + P9 y^3 + ...
491  Cx1 = P1 + P4 y + P8 y2 + ...
492  Cx2 = P3 + P7 y + ...
493  Cx3 = P6 + ...
494  ...
495 
496  Compute Cx0, Cx1...Cxn by solving 1-d polynomials in y in the usual way.
497  These values are cached and only recomputed for new values of Y or if the parameters change.
498 
499  Then compute f(x,y) by solving the 1-d polynomial in x in the usual way.
500  */
501  const int maxXCoeffInd = this->_order;
502 
503  if ((y != _oldY) || !this->_isCacheValid) {
504  // update _xCoeffs cache
505  // note: paramInd is decremented in both of the following loops
506  int paramInd = static_cast<int>(this->_params.size()) - 1;
507 
508  // initialize _xCoeffs to coeffs for pure y^n; e.g. for 3rd order:
509  // _xCoeffs[0] = _params[9], _xCoeffs[1] = _params[8], ... _xCoeffs[3] = _params[6]
510  for (int xCoeffInd = 0; xCoeffInd <= maxXCoeffInd; ++xCoeffInd, --paramInd) {
511  _xCoeffs[xCoeffInd] = this->_params[paramInd];
512  }
513 
514  // finish computing _xCoeffs
515  for (int xCoeffInd = 0, endXCoeffInd = maxXCoeffInd; paramInd >= 0; --paramInd) {
516  _xCoeffs[xCoeffInd] = (_xCoeffs[xCoeffInd] * y) + this->_params[paramInd];
517  ++xCoeffInd;
518  if (xCoeffInd >= endXCoeffInd) {
519  xCoeffInd = 0;
520  --endXCoeffInd;
521  }
522  }
523 
524  _oldY = y;
525  this->_isCacheValid = true;
526  }
527 
528  // use _xCoeffs to compute result
529  double retVal = _xCoeffs[maxXCoeffInd];
530  for (int xCoeffInd = maxXCoeffInd - 1; xCoeffInd >= 0; --xCoeffInd) {
531  retVal = (retVal * x) + _xCoeffs[xCoeffInd];
532  }
533  return static_cast<ReturnT>(retVal);
534  }
535 
540  std::vector<double> getDFuncDParameters(double x, double y) const override;
541 
542  std::string toString(std::string const& prefix) const override {
544  os << "PolynomialFunction2 [" << this->_order << "]: ";
545  os << Function2<ReturnT>::toString(prefix);
546  return os.str();
547  }
548 
549  bool isPersistable() const noexcept override { return true; }
550 
551 protected:
552  std::string getPersistenceName() const override;
553 
554  void write(afw::table::io::OutputArchiveHandle& handle) const override;
555 
556 private:
557  mutable double _oldY;
558  mutable std::vector<double> _xCoeffs;
559 
560 protected:
561  /* Default constructor: intended only for serialization */
562  explicit PolynomialFunction2() : BasePolynomialFunction2<ReturnT>(), _oldY(0), _xCoeffs(0) {}
563 };
564 
581 template <typename ReturnT>
582 class Chebyshev1Function1 : public Function1<ReturnT> {
583 public:
589  explicit Chebyshev1Function1(unsigned int order,
590  double minX = -1,
591  double maxX = 1)
592  : Function1<ReturnT>(order + 1) {
593  _initialize(minX, maxX);
594  }
595 
604  double minX = -1,
605  double maxX = 1)
606  : Function1<ReturnT>(params) {
607  if (params.size() < 1) {
609  "Chebyshev1Function1 called with empty vector");
610  }
611  _initialize(minX, maxX);
612  }
613 
618  ~Chebyshev1Function1() noexcept override = default;
619 
620  std::shared_ptr<Function1<ReturnT>> clone() const override {
621  return std::shared_ptr<Function1<ReturnT>>(new Chebyshev1Function1(this->_params, _minX, _maxX));
622  }
623 
627  double getMinX() const noexcept { return _minX; };
628 
632  double getMaxX() const noexcept { return _maxX; };
633 
637  unsigned int getOrder() const noexcept { return this->getNParameters() - 1; };
638 
639  bool isLinearCombination() const noexcept override { return true; };
640 
641  ReturnT operator()(double x) const override {
642  double xPrime = (x + _offset) * _scale;
643 
644  // Clenshaw function for solving the Chebyshev polynomial
645  // Non-recursive version from Kresimir Cosic
646  int const order = _order;
647  if (order == 0) {
648  return this->_params[0];
649  } else if (order == 1) {
650  return this->_params[0] + (this->_params[1] * xPrime);
651  }
652  double cshPrev = this->_params[order];
653  double csh = (2 * xPrime * this->_params[order]) + this->_params[order - 1];
654  for (int i = order - 2; i > 0; --i) {
655  double cshNext = (2 * xPrime * csh) + this->_params[i] - cshPrev;
656  cshPrev = csh;
657  csh = cshNext;
658  }
659  return (xPrime * csh) + this->_params[0] - cshPrev;
660  }
661 
662  std::string toString(std::string const& prefix) const override {
664  os << "Chebyshev1Function1 [" << _minX << ", " << _maxX << "]: ";
665  os << Function1<ReturnT>::toString(prefix);
666  return os.str();
667  }
668 
669 private:
670  double _minX;
671  double _maxX;
672  double _scale;
673  double _offset;
674  unsigned int _order;
675 
679  void _initialize(double minX, double maxX) {
680  _minX = minX;
681  _maxX = maxX;
682  _scale = 2.0 / (_maxX - _minX);
683  _offset = -(_minX + _maxX) * 0.5;
684  _order = this->getNParameters() - 1;
685  }
686 
687 protected:
688  /* Default constructor: intended only for serialization */
690  : Function1<ReturnT>(1), _minX(0.0), _maxX(0.0), _scale(1.0), _offset(0.0), _order(0) {}
691 };
692 
718 template <typename ReturnT>
720 public:
726  explicit Chebyshev1Function2(unsigned int order,
727  lsst::geom::Box2D const& xyRange = lsst::geom::Box2D(
728  lsst::geom::Point2D(-1.0, -1.0),
729  lsst::geom::Point2D(1.0, 1.0)))
730  : BasePolynomialFunction2<ReturnT>(order),
731  _oldYPrime(0),
732  _yCheby(this->_order + 1),
733  _xCoeffs(this->_order + 1) {
734  _initialize(xyRange);
735  }
736 
746  lsst::geom::Box2D const& xyRange = lsst::geom::Box2D(
747  lsst::geom::Point2D(-1.0, -1.0),
748  lsst::geom::Point2D(1.0, 1.0)))
749  : BasePolynomialFunction2<ReturnT>(params),
750  _oldYPrime(0),
751  _yCheby(this->_order + 1),
752  _xCoeffs(this->_order + 1) {
753  _initialize(xyRange);
754  }
755 
760  ~Chebyshev1Function2() noexcept override = default;
761 
762  std::shared_ptr<Function2<ReturnT>> clone() const override {
764  new Chebyshev1Function2(this->_params, this->getXYRange()));
765  }
766 
770  lsst::geom::Box2D getXYRange() const noexcept {
771  return lsst::geom::Box2D(lsst::geom::Point2D(_minX, _minY), lsst::geom::Point2D(_maxX, _maxY));
772  };
773 
779  virtual Chebyshev1Function2 truncate(int truncOrder
780  ) const {
781  if (truncOrder > this->_order) {
783  os << "truncated order=" << truncOrder << " must be <= original order=" << this->_order;
785  }
786  int truncNParams = this->nParametersFromOrder(truncOrder);
787  std::vector<double> truncParams(this->_params.begin(), this->_params.begin() + truncNParams);
788  return Chebyshev1Function2(truncParams, this->getXYRange());
789  }
790 
791  ReturnT operator()(double x, double y) const override {
792  /* Solve as follows:
793  - f(x,y) = Cy0 T0(y') + Cy1 T1(y') + Cy2 T2(y') + Cy3 T3(y') + ...
794  where:
795  Cy0 = P0 T0(x') + P1 T1(x') + P3 T2(x') + P6 T3(x') + ...
796  Cy1 = P2 T0(x') + P4 T1(x') + P7 T2(x') + ...
797  Cy2 = P5 T0(x') + P8 T1(x') + ...
798  Cy3 = P9 T0(x') + ...
799  ...
800 
801  First compute Tn(x') for each n
802  Then use that to compute Cy0, Cy1, ...Cyn
803  Then solve the y Chebyshev polynomial using the Clenshaw algorithm
804  */
805  double const xPrime = (x + _offsetX) * _scaleX;
806  double const yPrime = (y + _offsetY) * _scaleY;
807 
808  const int nParams = static_cast<int>(this->_params.size());
809  const int order = this->_order;
810 
811  if (order == 0) {
812  return this->_params[0]; // No caching required
813  }
814 
815  if ((yPrime != _oldYPrime) || !this->_isCacheValid) {
816  // update cached _yCheby and _xCoeffs
817  _yCheby[0] = 1.0;
818  _yCheby[1] = yPrime;
819  for (int chebyInd = 2; chebyInd <= order; chebyInd++) {
820  _yCheby[chebyInd] = (2 * yPrime * _yCheby[chebyInd - 1]) - _yCheby[chebyInd - 2];
821  }
822 
823  for (int coeffInd = 0; coeffInd <= order; coeffInd++) {
824  _xCoeffs[coeffInd] = 0;
825  }
826  for (int coeffInd = 0, endCoeffInd = 0, paramInd = 0; paramInd < nParams; paramInd++) {
827  _xCoeffs[coeffInd] += this->_params[paramInd] * _yCheby[endCoeffInd];
828  --coeffInd;
829  ++endCoeffInd;
830  if (coeffInd < 0) {
831  coeffInd = endCoeffInd;
832  endCoeffInd = 0;
833  }
834  }
835 
836  _oldYPrime = yPrime;
837  this->_isCacheValid = true;
838  }
839 
840  // Clenshaw function for solving the Chebyshev polynomial
841  // Non-recursive version from Kresimir Cosic
842  if (order == 1) {
843  return _xCoeffs[0] + (_xCoeffs[1] * xPrime);
844  }
845  double cshPrev = _xCoeffs[order];
846  double csh = (2 * xPrime * _xCoeffs[order]) + _xCoeffs[order - 1];
847  for (int i = order - 2; i > 0; --i) {
848  double cshNext = (2 * xPrime * csh) + _xCoeffs[i] - cshPrev;
849  cshPrev = csh;
850  csh = cshNext;
851  }
852  return (xPrime * csh) + _xCoeffs[0] - cshPrev;
853  }
854 
855  std::string toString(std::string const& prefix) const override {
857  os << "Chebyshev1Function2 [";
858  os << this->_order << ", " << this->getXYRange() << "]";
859  os << Function2<ReturnT>::toString(prefix);
860  return os.str();
861  }
862 
863  bool isPersistable() const noexcept override { return true; }
864 
865 protected:
866  std::string getPersistenceName() const override;
867 
868  void write(afw::table::io::OutputArchiveHandle& handle) const override;
869 
870 private:
871  mutable double _oldYPrime;
872  mutable std::vector<double> _yCheby;
873  mutable std::vector<double> _xCoeffs;
874  double _minX;
875  double _minY;
876  double _maxX;
877  double _maxY;
878  double _scaleX;
879  double _scaleY;
880  double _offsetX;
881  double _offsetY;
882 
886  void _initialize(lsst::geom::Box2D const& xyRange) {
887  _minX = xyRange.getMinX();
888  _minY = xyRange.getMinY();
889  _maxX = xyRange.getMaxX();
890  _maxY = xyRange.getMaxY();
891  _scaleX = 2.0 / (_maxX - _minX);
892  _scaleY = 2.0 / (_maxY - _minY);
893  _offsetX = -(_minX + _maxX) * 0.5;
894  _offsetY = -(_minY + _maxY) * 0.5;
895  }
896 
897 protected:
898  /* Default constructor: intended only for serialization */
900  : BasePolynomialFunction2<ReturnT>(),
901  _oldYPrime(0),
902  _yCheby(0),
903  _xCoeffs(0),
904  _minX(0.0),
905  _minY(0.0),
906  _maxX(0.0),
907  _maxY(0.0),
908  _scaleX(1.0),
909  _scaleY(1.0),
910  _offsetX(0.0),
911  _offsetY(0.0) {}
912 };
913 
926 template <typename ReturnT>
927 class LanczosFunction1 : public Function1<ReturnT> {
928 public:
932  explicit LanczosFunction1(unsigned int n,
933  double xOffset = 0.0)
934  : Function1<ReturnT>(1), _invN(1.0 / static_cast<double>(n)) {
935  this->_params[0] = xOffset;
936  }
937 
942  ~LanczosFunction1() noexcept override = default;
943 
944  std::shared_ptr<Function1<ReturnT>> clone() const override {
946  }
947 
948  ReturnT operator()(double x) const noexcept(IS_NOTHROW_INIT<ReturnT>) override {
949  double xArg1 = (x - this->_params[0]) * lsst::geom::PI;
950  double xArg2 = xArg1 * _invN;
951  if (std::fabs(xArg1) > 1.0e-5) {
952  return static_cast<ReturnT>(std::sin(xArg1) * std::sin(xArg2) / (xArg1 * xArg2));
953  } else {
954  return static_cast<ReturnT>(1);
955  }
956  }
957 
961  unsigned int getOrder() const noexcept { return static_cast<unsigned int>(0.5 + (1.0 / _invN)); };
962 
963  std::string toString(std::string const& prefix) const override {
965  os << "LanczosFunction1 [" << this->getOrder() << "]: ";
966  ;
967  os << Function1<ReturnT>::toString(prefix);
968  return os.str();
969  }
970 
971 private:
972  double _invN; // == 1/n
973 
974 protected:
975  /* Default constructor: intended only for serialization */
976  explicit LanczosFunction1() : Function1<ReturnT>(1), _invN(1.0) {}
977 };
978 
991 template <typename ReturnT>
992 class LanczosFunction2 : public Function2<ReturnT> {
993 public:
997  explicit LanczosFunction2(unsigned int n,
998  double xOffset = 0.0,
999  double yOffset = 0.0)
1000  : Function2<ReturnT>(2), _invN(1.0 / static_cast<double>(n)) {
1001  this->_params[0] = xOffset;
1002  this->_params[1] = yOffset;
1003  }
1004 
1009  ~LanczosFunction2() noexcept override = default;
1010 
1011  std::shared_ptr<Function2<ReturnT>> clone() const override {
1013  new LanczosFunction2(this->getOrder(), this->_params[0], this->_params[1]));
1014  }
1015 
1016  ReturnT operator()(double x, double y) const noexcept(IS_NOTHROW_INIT<ReturnT>) override {
1017  double xArg1 = (x - this->_params[0]) * lsst::geom::PI;
1018  double xArg2 = xArg1 * _invN;
1019  double xFunc = 1;
1020  if (std::fabs(xArg1) > 1.0e-5) {
1021  xFunc = std::sin(xArg1) * std::sin(xArg2) / (xArg1 * xArg2);
1022  }
1023  double yArg1 = (y - this->_params[1]) * lsst::geom::PI;
1024  double yArg2 = yArg1 * _invN;
1025  double yFunc = 1;
1026  if (std::fabs(yArg1) > 1.0e-5) {
1027  yFunc = std::sin(yArg1) * std::sin(yArg2) / (yArg1 * yArg2);
1028  }
1029  return static_cast<ReturnT>(xFunc * yFunc);
1030  }
1031 
1035  unsigned int getOrder() const noexcept { return static_cast<unsigned int>(0.5 + (1.0 / _invN)); };
1036 
1037  std::string toString(std::string const& prefix) const override {
1039  os << "LanczosFunction2 [" << this->getOrder() << "]: ";
1040  ;
1041  os << Function2<ReturnT>::toString(prefix);
1042  return os.str();
1043  }
1044 
1045 private:
1046  double _invN;
1047 
1048 protected:
1049  /* Default constructor: intended only for serialization */
1050  explicit LanczosFunction2() : Function2<ReturnT>(2), _invN(1.0) {}
1051 };
1052 } // namespace math
1053 } // namespace afw
1054 } // namespace lsst
1055 
1056 #endif // #ifndef LSST_AFW_MATH_FUNCTIONLIBRARY_H
double x
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
table::Key< double > angle
table::Key< double > ampl2
table::Key< double > sigma2
table::Key< double > sigma1
afw::table::Key< double > sigma
Definition: GaussianPsf.cc:49
std::ostream * os
Definition: Schema.cc:557
std::string prefix
Definition: SchemaMapper.cc:72
int y
Definition: SpanSet.cc:48
table::Key< int > b
T begin(T... args)
Base class for 2-dimensional polynomials of the form:
Definition: Function.h:326
static int nParametersFromOrder(int order)
Compute number of parameters from polynomial order.
Definition: Function.h:369
1-dimensional weighted sum of Chebyshev polynomials of the first kind.
double getMinX() const noexcept
Get minimum allowed x.
Chebyshev1Function1(std::vector< double > params, double minX=-1, double maxX=1)
Construct a Chebyshev polynomial with specified parameters and range.
Chebyshev1Function1 & operator=(Chebyshev1Function1 const &)=default
unsigned int getOrder() const noexcept
Get the polynomial order.
bool isLinearCombination() const noexcept override
Is the function a linear combination of its parameters?
Chebyshev1Function1(Chebyshev1Function1 const &)=default
Chebyshev1Function1(Chebyshev1Function1 &&)=default
ReturnT operator()(double x) const override
Chebyshev1Function1(unsigned int order, double minX=-1, double maxX=1)
Construct a Chebyshev polynomial of specified order and range.
double getMaxX() const noexcept
Get maximum allowed x.
~Chebyshev1Function1() noexcept override=default
std::shared_ptr< Function1< ReturnT > > clone() const override
Return a pointer to a deep copy of this function.
Chebyshev1Function1 & operator=(Chebyshev1Function1 &&)=default
std::string toString(std::string const &prefix) const override
Return a string representation of the function.
2-dimensional weighted sum of Chebyshev polynomials of the first kind.
std::string toString(std::string const &prefix) const override
Return a string representation of the function.
Chebyshev1Function2 & operator=(Chebyshev1Function2 const &)=default
ReturnT operator()(double x, double y) const override
bool isPersistable() const noexcept override
Return true if this particular object can be persisted using afw::table::io.
std::shared_ptr< Function2< ReturnT > > clone() const override
Return a pointer to a deep copy of this function.
Chebyshev1Function2(unsigned int order, lsst::geom::Box2D const &xyRange=lsst::geom::Box2D(lsst::geom::Point2D(-1.0, -1.0), lsst::geom::Point2D(1.0, 1.0)))
Construct a Chebyshev polynomial of specified order and range.
Chebyshev1Function2(Chebyshev1Function2 const &)=default
~Chebyshev1Function2() noexcept override=default
lsst::geom::Box2D getXYRange() const noexcept
Get x,y range.
Chebyshev1Function2 & operator=(Chebyshev1Function2 &&)=default
std::string getPersistenceName() const override
Return the unique name used to persist this object and look up its factory.
virtual Chebyshev1Function2 truncate(int truncOrder) const
Return a truncated copy of lower (or equal) order.
void write(afw::table::io::OutputArchiveHandle &handle) const override
Write the object to one or more catalogs.
Chebyshev1Function2(Chebyshev1Function2 &&)=default
Chebyshev1Function2(std::vector< double > params, lsst::geom::Box2D const &xyRange=lsst::geom::Box2D(lsst::geom::Point2D(-1.0, -1.0), lsst::geom::Point2D(1.0, 1.0)))
Construct a Chebyshev polynomial with specified parameters and range.
double Guassian (sum of two Gaussians)
void write(afw::table::io::OutputArchiveHandle &handle) const override
Write the object to one or more catalogs.
~DoubleGaussianFunction2() noexcept override=default
DoubleGaussianFunction2(double sigma1, double sigma2=0, double ampl2=0)
Construct a Gaussian function with specified x and y sigma.
DoubleGaussianFunction2 & operator=(DoubleGaussianFunction2 const &)=default
std::shared_ptr< Function2< ReturnT > > clone() const override
Return a pointer to a deep copy of this function.
DoubleGaussianFunction2(DoubleGaussianFunction2 &&)=default
DoubleGaussianFunction2(DoubleGaussianFunction2 const &)=default
std::string toString(std::string const &prefix) const override
Return a string representation of the function.
bool isPersistable() const noexcept override
Return true if this particular object can be persisted using afw::table::io.
ReturnT operator()(double x, double y) const noexcept(IS_NOTHROW_INIT< ReturnT >) override
DoubleGaussianFunction2 & operator=(DoubleGaussianFunction2 &&)=default
std::string getPersistenceName() const override
Return the unique name used to persist this object and look up its factory.
A Function taking one argument.
Definition: Function.h:202
A Function taking two arguments.
Definition: Function.h:259
std::vector< double > _params
Definition: Function.h:185
unsigned int getNParameters() const noexcept
Return the number of function parameters.
Definition: Function.h:112
~GaussianFunction1() noexcept override=default
GaussianFunction1(GaussianFunction1 &&)=default
std::string toString(std::string const &prefix) const override
Return a string representation of the function.
GaussianFunction1 & operator=(GaussianFunction1 &&)=default
GaussianFunction1(GaussianFunction1 const &)=default
GaussianFunction1(double sigma)
Construct a Gaussian function with specified sigma.
GaussianFunction1 & operator=(GaussianFunction1 const &)=default
std::shared_ptr< Function1< ReturnT > > clone() const override
Return a pointer to a deep copy of this function.
ReturnT operator()(double x) const noexcept(IS_NOTHROW_INIT< ReturnT >) override
GaussianFunction2(double sigma1, double sigma2, double angle=0.0)
Construct a 2-dimensional Gaussian function.
GaussianFunction2(GaussianFunction2 &&)=default
GaussianFunction2(GaussianFunction2 const &)=default
std::string toString(std::string const &prefix) const override
Return a string representation of the function.
~GaussianFunction2() noexcept override=default
void write(afw::table::io::OutputArchiveHandle &handle) const override
Write the object to one or more catalogs.
GaussianFunction2 & operator=(GaussianFunction2 const &)=default
GaussianFunction2 & operator=(GaussianFunction2 &&)=default
ReturnT operator()(double x, double y) const noexcept(IS_NOTHROW_INIT< ReturnT >) override
std::shared_ptr< Function2< ReturnT > > clone() const override
Return a pointer to a deep copy of this function.
std::string getPersistenceName() const override
Return the unique name used to persist this object and look up its factory.
bool isPersistable() const noexcept override
Return true if this particular object can be persisted using afw::table::io.
1-dimensional integer delta function.
~IntegerDeltaFunction1() noexcept override=default
IntegerDeltaFunction1(IntegerDeltaFunction1 const &)=default
ReturnT operator()(double x) const noexcept(IS_NOTHROW_INIT< ReturnT >) override
IntegerDeltaFunction1 & operator=(IntegerDeltaFunction1 const &)=default
IntegerDeltaFunction1(double xo)
Construct an integer delta function with specified xo.
IntegerDeltaFunction1 & operator=(IntegerDeltaFunction1 &&)=default
IntegerDeltaFunction1(IntegerDeltaFunction1 &&)=default
std::string toString(std::string const &prefix="") const override
Return a string representation of the function.
std::shared_ptr< Function1< ReturnT > > clone() const override
Return a pointer to a deep copy of this function.
2-dimensional integer delta function.
std::string toString(std::string const &prefix) const override
Return a string representation of the function.
IntegerDeltaFunction2(IntegerDeltaFunction2 const &)=default
IntegerDeltaFunction2 & operator=(IntegerDeltaFunction2 const &)=default
ReturnT operator()(double x, double y) const noexcept(IS_NOTHROW_INIT< ReturnT >) override
~IntegerDeltaFunction2() noexcept override=default
IntegerDeltaFunction2(double xo, double yo)
Construct an integer delta function with specified xo, yo.
IntegerDeltaFunction2 & operator=(IntegerDeltaFunction2 &&)=default
IntegerDeltaFunction2(IntegerDeltaFunction2 &&)=default
std::shared_ptr< Function2< ReturnT > > clone() const override
Return a pointer to a deep copy of this function.
1-dimensional Lanczos function
LanczosFunction1(unsigned int n, double xOffset=0.0)
Construct a Lanczos function of specified order and x,y offset.
~LanczosFunction1() noexcept override=default
LanczosFunction1 & operator=(LanczosFunction1 &&)=default
LanczosFunction1(LanczosFunction1 &&)=default
std::shared_ptr< Function1< ReturnT > > clone() const override
Return a pointer to a deep copy of this function.
LanczosFunction1(LanczosFunction1 const &)=default
ReturnT operator()(double x) const noexcept(IS_NOTHROW_INIT< ReturnT >) override
LanczosFunction1 & operator=(LanczosFunction1 const &)=default
unsigned int getOrder() const noexcept
Get the order of the Lanczos function.
std::string toString(std::string const &prefix) const override
Return a string representation of the function.
2-dimensional separable Lanczos function
std::shared_ptr< Function2< ReturnT > > clone() const override
Return a pointer to a deep copy of this function.
~LanczosFunction2() noexcept override=default
LanczosFunction2 & operator=(LanczosFunction2 &&)=default
LanczosFunction2(LanczosFunction2 &&)=default
LanczosFunction2(LanczosFunction2 const &)=default
LanczosFunction2(unsigned int n, double xOffset=0.0, double yOffset=0.0)
Construct a Lanczos function of specified order and x,y offset.
ReturnT operator()(double x, double y) const noexcept(IS_NOTHROW_INIT< ReturnT >) override
unsigned int getOrder() const noexcept
Get the order of Lanczos function.
LanczosFunction2 & operator=(LanczosFunction2 const &)=default
std::string toString(std::string const &prefix) const override
Return a string representation of the function.
1-dimensional polynomial function.
unsigned int getOrder() const noexcept
Get the polynomial order.
PolynomialFunction1 & operator=(PolynomialFunction1 &&)=default
std::shared_ptr< Function1< ReturnT > > clone() const override
Return a pointer to a deep copy of this function.
std::string toString(std::string const &prefix) const override
Return a string representation of the function.
bool isLinearCombination() const noexcept override
Is the function a linear combination of its parameters?
PolynomialFunction1(std::vector< double > params)
Construct a polynomial function with the specified parameters.
PolynomialFunction1 & operator=(PolynomialFunction1 const &)=default
PolynomialFunction1(PolynomialFunction1 &&)=default
PolynomialFunction1(unsigned int order)
Construct a polynomial function of the specified order.
PolynomialFunction1(PolynomialFunction1 const &)=default
ReturnT operator()(double x) const noexcept(IS_NOTHROW_INIT< ReturnT >) override
~PolynomialFunction1() noexcept override=default
2-dimensional polynomial function with cross terms
PolynomialFunction2(std::vector< double > params)
Construct a polynomial function with specified parameters.
PolynomialFunction2 & operator=(PolynomialFunction2 const &)=default
PolynomialFunction2(unsigned int order)
Construct a polynomial function of specified order.
PolynomialFunction2 & operator=(PolynomialFunction2 &&)=default
std::string toString(std::string const &prefix) const override
Return a string representation of the function.
PolynomialFunction2(PolynomialFunction2 const &)=default
std::shared_ptr< Function2< ReturnT > > clone() const override
Return a pointer to a deep copy of this function.
~PolynomialFunction2() noexcept override=default
std::vector< double > getDFuncDParameters(double x, double y) const override
Return the coefficients of the Function's parameters, evaluated at (x, y) I.e.
Definition: Function.cc:34
bool isPersistable() const noexcept override
Return true if this particular object can be persisted using afw::table::io.
std::string getPersistenceName() const override
Return the unique name used to persist this object and look up its factory.
ReturnT operator()(double x, double y) const noexcept(IS_NOTHROW_INIT< ReturnT >) override
PolynomialFunction2(PolynomialFunction2 &&)=default
void write(afw::table::io::OutputArchiveHandle &handle) const override
Write the object to one or more catalogs.
An object passed to Persistable::write to allow it to persist itself.
A floating-point coordinate rectangle geometry.
Definition: Box.h:413
double getMaxY() const noexcept
Definition: Box.h:519
double getMaxX() const noexcept
Definition: Box.h:518
double getMinY() const noexcept
Definition: Box.h:515
double getMinX() const noexcept
Definition: Box.h:514
Reports invalid arguments.
Definition: Runtime.h:66
T cos(T... args)
T exp(T... args)
T fabs(T... args)
constexpr double PI
The ratio of a circle's circumference to diameter.
Definition: Angle.h:39
constexpr double TWOPI
Definition: Angle.h:40
A base class for image defects.
STL namespace.
T sin(T... args)
T size(T... args)
table::Key< int > order