LSSTApplications  17.0+124,17.0+14,17.0+73,18.0.0+37,18.0.0+80,18.0.0-4-g68ffd23+4,18.1.0-1-g0001055+12,18.1.0-1-g03d53ef+5,18.1.0-1-g1349e88+55,18.1.0-1-g2505f39+44,18.1.0-1-g5315e5e+4,18.1.0-1-g5e4b7ea+14,18.1.0-1-g7e8fceb+4,18.1.0-1-g85f8cd4+48,18.1.0-1-g8ff0b9f+4,18.1.0-1-ga2c679d+1,18.1.0-1-gd55f500+35,18.1.0-10-gb58edde+2,18.1.0-11-g0997b02+4,18.1.0-13-gfe4edf0b+12,18.1.0-14-g259bd21+21,18.1.0-19-gdb69f3f+2,18.1.0-2-g5f9922c+24,18.1.0-2-gd3b74e5+11,18.1.0-2-gfbf3545+32,18.1.0-26-g728bddb4+5,18.1.0-27-g6ff7ca9+2,18.1.0-3-g52aa583+25,18.1.0-3-g8ea57af+9,18.1.0-3-gb69f684+42,18.1.0-3-gfcaddf3+6,18.1.0-32-gd8786685a,18.1.0-4-gf3f9b77+6,18.1.0-5-g1dd662b+2,18.1.0-5-g6dbcb01+41,18.1.0-6-gae77429+3,18.1.0-7-g9d75d83+9,18.1.0-7-gae09a6d+30,18.1.0-9-gc381ef5+4,w.2019.45
LSSTDataManagementBasePackage
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 
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 
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  }
156  GaussianFunction1(GaussianFunction1 const&) = default;
158  GaussianFunction1& operator=(GaussianFunction1 const&) = default;
160  ~GaussianFunction1() noexcept override = default;
161 
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 
216  GaussianFunction2(GaussianFunction2 const&) = default;
218  GaussianFunction2& operator=(GaussianFunction2 const&) = default;
220  ~GaussianFunction2() noexcept override = default;
221 
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 
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 
396  PolynomialFunction1(PolynomialFunction1 const&) = default;
400  ~PolynomialFunction1() noexcept override = default;
401 
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)
473  : BasePolynomialFunction2<ReturnT>(params), _oldY(0), _xCoeffs(this->_order + 1) {}
475 
476  PolynomialFunction2(PolynomialFunction2 const&) = default;
480  ~PolynomialFunction2() noexcept override = default;
481 
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 
614  Chebyshev1Function1(Chebyshev1Function1 const&) = default;
618  ~Chebyshev1Function1() noexcept override = default;
619 
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 
745  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 
756  Chebyshev1Function2(Chebyshev1Function2 const&) = default;
760  ~Chebyshev1Function2() noexcept override = default;
761 
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 
938  LanczosFunction1(LanczosFunction1 const&) = default;
939  LanczosFunction1(LanczosFunction1&&) = default;
940  LanczosFunction1& operator=(LanczosFunction1 const&) = default;
942  ~LanczosFunction1() noexcept override = default;
943 
945  return std::shared_ptr<Function1<ReturnT>>(new LanczosFunction1(this->getOrder(), this->_params[0]));
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 
1005  LanczosFunction2(LanczosFunction2 const&) = default;
1006  LanczosFunction2(LanczosFunction2&&) = default;
1007  LanczosFunction2& operator=(LanczosFunction2 const&) = default;
1009  ~LanczosFunction2() noexcept override = default;
1010 
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 getMinY() const noexcept
Definition: Box.h:515
ReturnT operator()(double x, double y) const noexcept(IS_NOTHROW_INIT< ReturnT >) override
2-dimensional weighted sum of Chebyshev polynomials of the first kind.
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 override
ReturnT operator()(double x, double y) const noexcept(IS_NOTHROW_INIT< ReturnT >) override
unsigned int getNParameters() const noexcept
Return the number of function parameters.
Definition: Function.h:112
std::shared_ptr< Function2< ReturnT > > clone() const override
Return a pointer to a deep copy of this function.
PolynomialFunction1(unsigned int order)
Construct a polynomial function of the specified order.
ReturnT operator()(double x, double y) const noexcept(IS_NOTHROW_INIT< ReturnT >) override
PolynomialFunction2(unsigned int order)
Construct a polynomial function of specified order.
std::shared_ptr< Function1< ReturnT > > clone() const override
Return a pointer to a deep copy of this function.
1-dimensional Lanczos function
double Guassian (sum of two Gaussians)
double getMaxX() const noexcept
Get maximum allowed x.
unsigned int getOrder() const noexcept
Get the order of the Lanczos function.
A floating-point coordinate rectangle geometry.
Definition: Box.h:413
An object passed to Persistable::write to allow it to persist itself.
double constexpr TWOPI
Definition: Angle.h:40
std::shared_ptr< Function2< ReturnT > > clone() const override
Return a pointer to a deep copy of this function.
GaussianFunction1(double sigma)
Construct a Gaussian function with specified sigma.
std::shared_ptr< Function2< ReturnT > > clone() const override
Return a pointer to a deep copy of this function.
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.
double getMinX() const noexcept
Definition: Box.h:514
T exp(T... args)
table::Key< int > b
bool isPersistable() const noexcept override
Return true if this particular object can be persisted using afw::table::io.
2-dimensional polynomial function with cross terms
PolynomialFunction1(std::vector< double > params)
Construct a polynomial function with the specified parameters.
1-dimensional weighted sum of Chebyshev polynomials of the first kind.
int y
Definition: SpanSet.cc:49
double getMaxX() const noexcept
Definition: Box.h:518
STL namespace.
Chebyshev1Function1(unsigned int order, double minX=-1, double maxX=1)
Construct a Chebyshev polynomial of specified order and range.
2-dimensional integer delta function.
virtual Chebyshev1Function2 truncate(int truncOrder) const
Return a truncated copy of lower (or equal) order.
IntegerDeltaFunction1(double xo)
Construct an integer delta function with specified xo.
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.
table::Key< double > angle
ReturnT operator()(double x, double y) const noexcept(IS_NOTHROW_INIT< ReturnT >) override
bool isLinearCombination() const noexcept override
Is the function a linear combination of its parameters?
std::shared_ptr< Function1< ReturnT > > clone() const override
Return a pointer to a deep copy of this function.
table::Key< double > sigma2
LanczosFunction2(unsigned int n, double xOffset=0.0, double yOffset=0.0)
Construct a Lanczos function of specified order and x,y offset.
IntegerDeltaFunction1 & operator=(IntegerDeltaFunction1 const &)=default
std::shared_ptr< Function2< ReturnT > > clone() const override
Return a pointer to a deep copy of this function.
std::shared_ptr< Function2< ReturnT > > clone() const override
Return a pointer to a deep copy of this function.
2-dimensional separable Lanczos function
A Function taking two arguments.
Definition: Function.h:259
bool isLinearCombination() const noexcept override
Is the function a linear combination of its parameters?
std::string toString(std::string const &prefix) const override
Return a string representation of the function.
std::string toString(std::string const &prefix) const override
Return a string representation of the function.
ReturnT operator()(double x) const override
virtual void write(OutputArchiveHandle &handle) const
Write the object to one or more catalogs.
Definition: Persistable.cc:38
Base class for 2-dimensional polynomials of the form:
Definition: Function.h:326
afw::table::Key< double > sigma
Definition: GaussianPsf.cc:50
STL class.
T sin(T... args)
double getMaxY() const noexcept
Definition: Box.h:519
std::vector< double > _params
Definition: Function.h:185
A base class for image defects.
IntegerDeltaFunction2(double xo, double yo)
Construct an integer delta function with specified xo, yo.
1-dimensional integer delta function.
unsigned int getOrder() const noexcept
Get the order of Lanczos function.
A Function taking one argument.
Definition: Function.h:202
ReturnT operator()(double x) const noexcept(IS_NOTHROW_INIT< ReturnT >) override
T str(T... args)
PolynomialFunction2(std::vector< double > params)
Construct a polynomial function with specified parameters.
T cos(T... args)
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.
ReturnT operator()(double x, double y) const noexcept(IS_NOTHROW_INIT< ReturnT >) override
T fabs(T... args)
ReturnT operator()(double x) const noexcept(IS_NOTHROW_INIT< ReturnT >) override
double x
std::shared_ptr< Function2< 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
T size(T... args)
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
std::string toString(std::string const &prefix) const override
Return a string representation of the function.
DoubleGaussianFunction2(double sigma1, double sigma2=0, double ampl2=0)
Construct a Gaussian function with specified x and y sigma.
std::string toString(std::string const &prefix) const override
Return a string representation of the function.
table::Key< double > ampl2
T begin(T... args)
lsst::geom::Box2D getXYRange() const noexcept
Get x,y range.
std::string toString(std::string const &prefix) const override
Return a string representation of the function.
unsigned int getOrder() const noexcept
Get the polynomial order.
~IntegerDeltaFunction1() noexcept override=default
ReturnT operator()(double x) const noexcept(IS_NOTHROW_INIT< ReturnT >) override
std::string toString(std::string const &prefix) const override
Return a string representation of the 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.
Reports invalid arguments.
Definition: Runtime.h:66
std::string prefix
Definition: SchemaMapper.cc:79
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.
virtual std::string getPersistenceName() const
Return the unique name used to persist this object and look up its factory.
Definition: Persistable.cc:34
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.
GaussianFunction2(double sigma1, double sigma2, double angle=0.0)
Construct a 2-dimensional Gaussian function.
bool isPersistable() const noexcept override
Return true if this particular object can be persisted using afw::table::io.
std::string toString(std::string const &prefix="") const override
Return a string representation of the function.
std::ostream * os
Definition: Schema.cc:746
bool isPersistable() const noexcept override
Return true if this particular object can be persisted using afw::table::io.
double constexpr PI
The ratio of a circle&#39;s circumference to diameter.
Definition: Angle.h:39
unsigned int getOrder() const noexcept
Get the polynomial order.
std::string toString(std::string const &prefix) const override
Return a string representation of the function.
table::Key< double > sigma1
1-dimensional polynomial function.
LanczosFunction1(unsigned int n, double xOffset=0.0)
Construct a Lanczos function of specified order and x,y offset.