LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
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"
35#include "lsst/geom/Angle.h"
36
37namespace lsst {
38namespace afw {
39namespace math {
40
50template <typename ReturnT>
51class IntegerDeltaFunction1 : public Function1<ReturnT> {
52public:
56 explicit IntegerDeltaFunction1(double xo) : Function1<ReturnT>(0), _xo(xo) {}
57
62 ~IntegerDeltaFunction1() noexcept override = default;
63
64 std::shared_ptr<Function1<ReturnT>> clone() const override {
66 }
67
68 ReturnT operator()(double x) const noexcept(IS_NOTHROW_INIT<ReturnT>) override {
69 return static_cast<ReturnT>(x == _xo);
70 }
71
72 std::string toString(std::string const& prefix = "") const override {
74 os << "IntegerDeltaFunction1 [" << _xo << "]: ";
75 os << Function1<ReturnT>::toString(prefix);
76 return os.str();
77 };
78
79private:
80 double _xo;
81
82protected:
83 /* Default constructor: intended only for serialization */
84 explicit IntegerDeltaFunction1() : Function1<ReturnT>(0), _xo(0.0) {}
85};
86
96template <typename ReturnT>
97class IntegerDeltaFunction2 : public Function2<ReturnT> {
98public:
102 explicit IntegerDeltaFunction2(double xo, double yo) : Function2<ReturnT>(0), _xo(xo), _yo(yo) {}
103
108 ~IntegerDeltaFunction2() noexcept override = default;
109
110 std::shared_ptr<Function2<ReturnT>> clone() const override {
112 }
113
114 ReturnT operator()(double x, double y) const noexcept(IS_NOTHROW_INIT<ReturnT>) override {
115 return static_cast<ReturnT>((x == _xo) && (y == _yo));
116 }
117
118 std::string toString(std::string const& prefix) const override {
120 os << "IntegerDeltaFunction2 [" << _xo << ", " << _yo << "]: ";
121 os << Function2<ReturnT>::toString(prefix);
122 return os.str();
123 }
124
125private:
126 double _xo;
127 double _yo;
128
129protected:
130 /* Default constructor: intended only for serialization */
131 explicit IntegerDeltaFunction2() : Function2<ReturnT>(0), _xo(0.0), _yo(0.0) {}
132};
133
144template <typename ReturnT>
145class GaussianFunction1 : public Function1<ReturnT> {
146public:
150 explicit GaussianFunction1(double sigma)
151 : Function1<ReturnT>(1), _multFac(1.0 / std::sqrt(lsst::geom::TWOPI)) {
152 this->_params[0] = sigma;
153 }
158 ~GaussianFunction1() noexcept override = default;
159
160 std::shared_ptr<Function1<ReturnT>> clone() const override {
162 }
163
164 ReturnT operator()(double x) const noexcept(IS_NOTHROW_INIT<ReturnT>) override {
165 return static_cast<ReturnT>((_multFac / this->_params[0]) *
166 std::exp(-(x * x) / (2.0 * this->_params[0] * this->_params[0])));
167 }
168
169 std::string toString(std::string const& prefix) const override {
171 os << "GaussianFunction1 [" << _multFac << "]: ";
172 os << Function1<ReturnT>::toString(prefix);
173 return os.str();
174 }
175
176private:
177 const double _multFac;
178
179protected:
180 /* Default constructor: intended only for serialization */
181 explicit GaussianFunction1() : Function1<ReturnT>(1), _multFac(1.0 / std::sqrt(lsst::geom::TWOPI)) {}
182};
183
198template <typename ReturnT>
199class GaussianFunction2 : public Function2<ReturnT> {
200public:
204 explicit GaussianFunction2(double sigma1,
205 double sigma2,
206 double angle = 0.0)
207 : Function2<ReturnT>(3), _multFac(1.0 / (lsst::geom::TWOPI)) {
208 this->_params[0] = sigma1;
209 this->_params[1] = sigma2;
210 this->_params[2] = angle;
211 _updateCache();
212 }
213
218 ~GaussianFunction2() noexcept override = default;
219
220 std::shared_ptr<Function2<ReturnT>> clone() const override {
222 new GaussianFunction2(this->_params[0], this->_params[1], this->_params[2]));
223 }
224
225 ReturnT operator()(double x, double y) const noexcept(IS_NOTHROW_INIT<ReturnT>) override {
226 if (_angle != this->_params[2]) {
227 _updateCache();
228 }
229 double pos1 = (_cosAngle * x) + (_sinAngle * y);
230 double pos2 = (-_sinAngle * x) + (_cosAngle * y);
231 return static_cast<ReturnT>((_multFac / (this->_params[0] * this->_params[1])) *
232 std::exp(-((pos1 * pos1) / (2.0 * this->_params[0] * this->_params[0])) -
233 ((pos2 * pos2) / (2.0 * this->_params[1] * this->_params[1]))));
234 }
235
236 std::string toString(std::string const& prefix) const override {
238 os << "GaussianFunction2: ";
239 os << Function2<ReturnT>::toString(prefix);
240 return os.str();
241 }
242
243 bool isPersistable() const noexcept override { return true; }
244
245protected:
246 std::string getPersistenceName() const override;
247
248 void write(afw::table::io::OutputArchiveHandle& handle) const override;
249
250private:
268 void _updateCache() const noexcept {
269 _angle = this->_params[2];
270 _sinAngle = std::sin(_angle);
271 _cosAngle = std::cos(_angle);
272 }
273 const double _multFac;
274 mutable double _angle;
275 mutable double _sinAngle;
276 mutable double _cosAngle;
277
278protected:
279 /* Default constructor: intended only for serialization */
281 : Function2<ReturnT>(3),
282 _multFac(1.0 / (lsst::geom::TWOPI)),
283 _angle(0.0),
284 _sinAngle(0.0),
285 _cosAngle(1.0) {}
286};
287
302template <typename ReturnT>
303class DoubleGaussianFunction2 : public Function2<ReturnT> {
304public:
309 double sigma1,
310 double sigma2 = 0,
311 double ampl2 = 0)
312 : Function2<ReturnT>(3), _multFac(1.0 / (lsst::geom::TWOPI)) {
313 this->_params[0] = sigma1;
314 this->_params[1] = sigma2;
315 this->_params[2] = ampl2;
316 }
317
322 ~DoubleGaussianFunction2() noexcept override = default;
323
324 std::shared_ptr<Function2<ReturnT>> clone() const override {
326 new DoubleGaussianFunction2(this->_params[0], this->_params[1], this->_params[2]));
327 }
328
329 ReturnT operator()(double x, double y) const noexcept(IS_NOTHROW_INIT<ReturnT>) override {
330 double radSq = (x * x) + (y * y);
331 double sigma1Sq = this->_params[0] * this->_params[0];
332 double sigma2Sq = this->_params[1] * this->_params[1];
333 double b = this->_params[2];
334 return static_cast<ReturnT>(
335 (_multFac / (sigma1Sq + (b * sigma2Sq))) *
336 (std::exp(-radSq / (2.0 * sigma1Sq)) + (b * std::exp(-radSq / (2.0 * sigma2Sq)))));
337 }
338
339 std::string toString(std::string const& prefix) const override {
341 os << "DoubleGaussianFunction2 [" << _multFac << "]: ";
342 os << Function2<ReturnT>::toString(prefix);
343 return os.str();
344 }
345
346 bool isPersistable() const noexcept override { return true; }
347
348protected:
349 std::string getPersistenceName() const override;
350
351 void write(afw::table::io::OutputArchiveHandle& handle) const override;
352
353private:
354 const double _multFac;
355
356protected:
357 /* Default constructor: intended only for serialization */
358 explicit DoubleGaussianFunction2() : Function2<ReturnT>(3), _multFac(1.0 / (lsst::geom::TWOPI)) {}
359};
360
368template <typename ReturnT>
369class PolynomialFunction1 : public Function1<ReturnT> {
370public:
376 explicit PolynomialFunction1(unsigned int order)
377 : Function1<ReturnT>(order + 1) {}
378
387 : Function1<ReturnT>(params) {
388 if (params.size() < 1) {
390 "PolynomialFunction1 called with empty vector");
391 }
392 }
393
398 ~PolynomialFunction1() noexcept override = default;
399
400 std::shared_ptr<Function1<ReturnT>> clone() const override {
402 }
403
404 bool isLinearCombination() const noexcept override { return true; };
405
406 ReturnT operator()(double x) const noexcept(IS_NOTHROW_INIT<ReturnT>) override {
407 int const order = static_cast<int>(this->_params.size()) - 1;
408 double retVal = this->_params[order];
409 for (int ii = order - 1; ii >= 0; --ii) {
410 retVal = (retVal * x) + this->_params[ii];
411 }
412 return static_cast<ReturnT>(retVal);
413 }
414
418 unsigned int getOrder() const noexcept { return this->getNParameters() - 1; };
419
420 std::string toString(std::string const& prefix) const override {
422 os << "PolynomialFunction1 []: ";
423 os << Function1<ReturnT>::toString(prefix);
424 return os.str();
425 }
426
427protected:
428 /* Default constructor: intended only for serialization */
429 explicit PolynomialFunction1() : Function1<ReturnT>(1) {}
430};
431
446template <typename ReturnT>
448public:
456 explicit PolynomialFunction2(unsigned int order)
457 : BasePolynomialFunction2<ReturnT>(order), _oldY(0), _xCoeffs(this->_order + 1) {}
458
470 std::vector<double> params)
472 : BasePolynomialFunction2<ReturnT>(params), _oldY(0), _xCoeffs(this->_order + 1) {}
473
478 ~PolynomialFunction2() noexcept override = default;
479
480 std::shared_ptr<Function2<ReturnT>> clone() const override {
482 }
483
484 ReturnT operator()(double x, double y) const noexcept(IS_NOTHROW_INIT<ReturnT>) override {
485 /* Solve as follows:
486 - f(x,y) = Cx0 + Cx1 x + Cx2 x^2 + Cx3 x^3 + ...
487 where:
488 Cx0 = P0 + P2 y + P5 y^2 + P9 y^3 + ...
489 Cx1 = P1 + P4 y + P8 y2 + ...
490 Cx2 = P3 + P7 y + ...
491 Cx3 = P6 + ...
492 ...
493
494 Compute Cx0, Cx1...Cxn by solving 1-d polynomials in y in the usual way.
495 These values are cached and only recomputed for new values of Y or if the parameters change.
496
497 Then compute f(x,y) by solving the 1-d polynomial in x in the usual way.
498 */
499 const int maxXCoeffInd = this->_order;
500
501 if ((y != _oldY) || !this->_isCacheValid) {
502 // update _xCoeffs cache
503 // note: paramInd is decremented in both of the following loops
504 int paramInd = static_cast<int>(this->_params.size()) - 1;
505
506 // initialize _xCoeffs to coeffs for pure y^n; e.g. for 3rd order:
507 // _xCoeffs[0] = _params[9], _xCoeffs[1] = _params[8], ... _xCoeffs[3] = _params[6]
508 for (int xCoeffInd = 0; xCoeffInd <= maxXCoeffInd; ++xCoeffInd, --paramInd) {
509 _xCoeffs[xCoeffInd] = this->_params[paramInd];
510 }
511
512 // finish computing _xCoeffs
513 for (int xCoeffInd = 0, endXCoeffInd = maxXCoeffInd; paramInd >= 0; --paramInd) {
514 _xCoeffs[xCoeffInd] = (_xCoeffs[xCoeffInd] * y) + this->_params[paramInd];
515 ++xCoeffInd;
516 if (xCoeffInd >= endXCoeffInd) {
517 xCoeffInd = 0;
518 --endXCoeffInd;
519 }
520 }
521
522 _oldY = y;
523 this->_isCacheValid = true;
524 }
525
526 // use _xCoeffs to compute result
527 double retVal = _xCoeffs[maxXCoeffInd];
528 for (int xCoeffInd = maxXCoeffInd - 1; xCoeffInd >= 0; --xCoeffInd) {
529 retVal = (retVal * x) + _xCoeffs[xCoeffInd];
530 }
531 return static_cast<ReturnT>(retVal);
532 }
533
538 std::vector<double> getDFuncDParameters(double x, double y) const override;
539
540 std::string toString(std::string const& prefix) const override {
542 os << "PolynomialFunction2 [" << this->_order << "]: ";
543 os << Function2<ReturnT>::toString(prefix);
544 return os.str();
545 }
546
547 bool isPersistable() const noexcept override { return true; }
548
549protected:
550 std::string getPersistenceName() const override;
551
552 void write(afw::table::io::OutputArchiveHandle& handle) const override;
553
554private:
555 mutable double _oldY;
556 mutable std::vector<double> _xCoeffs;
557
558protected:
559 /* Default constructor: intended only for serialization */
560 explicit PolynomialFunction2() : BasePolynomialFunction2<ReturnT>(), _oldY(0), _xCoeffs(0) {}
561};
562
579template <typename ReturnT>
580class Chebyshev1Function1 : public Function1<ReturnT> {
581public:
587 explicit Chebyshev1Function1(unsigned int order,
588 double minX = -1,
589 double maxX = 1)
590 : Function1<ReturnT>(order + 1) {
591 _initialize(minX, maxX);
592 }
593
602 double minX = -1,
603 double maxX = 1)
604 : Function1<ReturnT>(params) {
605 if (params.size() < 1) {
607 "Chebyshev1Function1 called with empty vector");
608 }
609 _initialize(minX, maxX);
610 }
611
616 ~Chebyshev1Function1() noexcept override = default;
617
618 std::shared_ptr<Function1<ReturnT>> clone() const override {
619 return std::shared_ptr<Function1<ReturnT>>(new Chebyshev1Function1(this->_params, _minX, _maxX));
620 }
621
625 double getMinX() const noexcept { return _minX; };
626
630 double getMaxX() const noexcept { return _maxX; };
631
635 unsigned int getOrder() const noexcept { return this->getNParameters() - 1; };
636
637 bool isLinearCombination() const noexcept override { return true; };
638
639 ReturnT operator()(double x) const override {
640 double xPrime = (x + _offset) * _scale;
641
642 // Clenshaw function for solving the Chebyshev polynomial
643 // Non-recursive version from Kresimir Cosic
644 int const order = _order;
645 if (order == 0) {
646 return this->_params[0];
647 } else if (order == 1) {
648 return this->_params[0] + (this->_params[1] * xPrime);
649 }
650 double cshPrev = this->_params[order];
651 double csh = (2 * xPrime * this->_params[order]) + this->_params[order - 1];
652 for (int i = order - 2; i > 0; --i) {
653 double cshNext = (2 * xPrime * csh) + this->_params[i] - cshPrev;
654 cshPrev = csh;
655 csh = cshNext;
656 }
657 return (xPrime * csh) + this->_params[0] - cshPrev;
658 }
659
660 std::string toString(std::string const& prefix) const override {
662 os << "Chebyshev1Function1 [" << _minX << ", " << _maxX << "]: ";
663 os << Function1<ReturnT>::toString(prefix);
664 return os.str();
665 }
666
667private:
668 double _minX;
669 double _maxX;
670 double _scale;
671 double _offset;
672 unsigned int _order;
673
677 void _initialize(double minX, double maxX) {
678 _minX = minX;
679 _maxX = maxX;
680 _scale = 2.0 / (_maxX - _minX);
681 _offset = -(_minX + _maxX) * 0.5;
682 _order = this->getNParameters() - 1;
683 }
684
685protected:
686 /* Default constructor: intended only for serialization */
688 : Function1<ReturnT>(1), _minX(0.0), _maxX(0.0), _scale(1.0), _offset(0.0), _order(0) {}
689};
690
716template <typename ReturnT>
718public:
724 explicit Chebyshev1Function2(unsigned int order,
725 lsst::geom::Box2D const& xyRange = lsst::geom::Box2D(
726 lsst::geom::Point2D(-1.0, -1.0),
727 lsst::geom::Point2D(1.0, 1.0)))
728 : BasePolynomialFunction2<ReturnT>(order),
729 _oldYPrime(0),
730 _yCheby(this->_order + 1),
731 _xCoeffs(this->_order + 1) {
732 _initialize(xyRange);
733 }
734
744 lsst::geom::Box2D const& xyRange = lsst::geom::Box2D(
745 lsst::geom::Point2D(-1.0, -1.0),
746 lsst::geom::Point2D(1.0, 1.0)))
747 : BasePolynomialFunction2<ReturnT>(params),
748 _oldYPrime(0),
749 _yCheby(this->_order + 1),
750 _xCoeffs(this->_order + 1) {
751 _initialize(xyRange);
752 }
753
758 ~Chebyshev1Function2() noexcept override = default;
759
760 std::shared_ptr<Function2<ReturnT>> clone() const override {
762 new Chebyshev1Function2(this->_params, this->getXYRange()));
763 }
764
768 lsst::geom::Box2D getXYRange() const noexcept {
769 return lsst::geom::Box2D(lsst::geom::Point2D(_minX, _minY), lsst::geom::Point2D(_maxX, _maxY));
770 };
771
777 virtual Chebyshev1Function2 truncate(int truncOrder
778 ) const {
779 if (truncOrder > this->_order) {
781 os << "truncated order=" << truncOrder << " must be <= original order=" << this->_order;
783 }
784 int truncNParams = this->nParametersFromOrder(truncOrder);
785 std::vector<double> truncParams(this->_params.begin(), this->_params.begin() + truncNParams);
786 return Chebyshev1Function2(truncParams, this->getXYRange());
787 }
788
789 ReturnT operator()(double x, double y) const override {
790 /* Solve as follows:
791 - f(x,y) = Cy0 T0(y') + Cy1 T1(y') + Cy2 T2(y') + Cy3 T3(y') + ...
792 where:
793 Cy0 = P0 T0(x') + P1 T1(x') + P3 T2(x') + P6 T3(x') + ...
794 Cy1 = P2 T0(x') + P4 T1(x') + P7 T2(x') + ...
795 Cy2 = P5 T0(x') + P8 T1(x') + ...
796 Cy3 = P9 T0(x') + ...
797 ...
798
799 First compute Tn(x') for each n
800 Then use that to compute Cy0, Cy1, ...Cyn
801 Then solve the y Chebyshev polynomial using the Clenshaw algorithm
802 */
803 double const xPrime = (x + _offsetX) * _scaleX;
804 double const yPrime = (y + _offsetY) * _scaleY;
805
806 const int nParams = static_cast<int>(this->_params.size());
807 const int order = this->_order;
808
809 if (order == 0) {
810 return this->_params[0]; // No caching required
811 }
812
813 if ((yPrime != _oldYPrime) || !this->_isCacheValid) {
814 // update cached _yCheby and _xCoeffs
815 _yCheby[0] = 1.0;
816 _yCheby[1] = yPrime;
817 for (int chebyInd = 2; chebyInd <= order; chebyInd++) {
818 _yCheby[chebyInd] = (2 * yPrime * _yCheby[chebyInd - 1]) - _yCheby[chebyInd - 2];
819 }
820
821 for (int coeffInd = 0; coeffInd <= order; coeffInd++) {
822 _xCoeffs[coeffInd] = 0;
823 }
824 for (int coeffInd = 0, endCoeffInd = 0, paramInd = 0; paramInd < nParams; paramInd++) {
825 _xCoeffs[coeffInd] += this->_params[paramInd] * _yCheby[endCoeffInd];
826 --coeffInd;
827 ++endCoeffInd;
828 if (coeffInd < 0) {
829 coeffInd = endCoeffInd;
830 endCoeffInd = 0;
831 }
832 }
833
834 _oldYPrime = yPrime;
835 this->_isCacheValid = true;
836 }
837
838 // Clenshaw function for solving the Chebyshev polynomial
839 // Non-recursive version from Kresimir Cosic
840 if (order == 1) {
841 return _xCoeffs[0] + (_xCoeffs[1] * xPrime);
842 }
843 double cshPrev = _xCoeffs[order];
844 double csh = (2 * xPrime * _xCoeffs[order]) + _xCoeffs[order - 1];
845 for (int i = order - 2; i > 0; --i) {
846 double cshNext = (2 * xPrime * csh) + _xCoeffs[i] - cshPrev;
847 cshPrev = csh;
848 csh = cshNext;
849 }
850 return (xPrime * csh) + _xCoeffs[0] - cshPrev;
851 }
852
853 std::string toString(std::string const& prefix) const override {
855 os << "Chebyshev1Function2 [";
856 os << this->_order << ", " << this->getXYRange() << "]";
857 os << Function2<ReturnT>::toString(prefix);
858 return os.str();
859 }
860
861 bool isPersistable() const noexcept override { return true; }
862
863protected:
864 std::string getPersistenceName() const override;
865
866 void write(afw::table::io::OutputArchiveHandle& handle) const override;
867
868private:
869 mutable double _oldYPrime;
870 mutable std::vector<double> _yCheby;
871 mutable std::vector<double> _xCoeffs;
872 double _minX;
873 double _minY;
874 double _maxX;
875 double _maxY;
876 double _scaleX;
877 double _scaleY;
878 double _offsetX;
879 double _offsetY;
880
884 void _initialize(lsst::geom::Box2D const& xyRange) {
885 _minX = xyRange.getMinX();
886 _minY = xyRange.getMinY();
887 _maxX = xyRange.getMaxX();
888 _maxY = xyRange.getMaxY();
889 _scaleX = 2.0 / (_maxX - _minX);
890 _scaleY = 2.0 / (_maxY - _minY);
891 _offsetX = -(_minX + _maxX) * 0.5;
892 _offsetY = -(_minY + _maxY) * 0.5;
893 }
894
895protected:
896 /* Default constructor: intended only for serialization */
898 : BasePolynomialFunction2<ReturnT>(),
899 _oldYPrime(0),
900 _yCheby(0),
901 _xCoeffs(0),
902 _minX(0.0),
903 _minY(0.0),
904 _maxX(0.0),
905 _maxY(0.0),
906 _scaleX(1.0),
907 _scaleY(1.0),
908 _offsetX(0.0),
909 _offsetY(0.0) {}
910};
911
924template <typename ReturnT>
925class LanczosFunction1 : public Function1<ReturnT> {
926public:
930 explicit LanczosFunction1(unsigned int n,
931 double xOffset = 0.0)
932 : Function1<ReturnT>(1), _invN(1.0 / static_cast<double>(n)) {
933 this->_params[0] = xOffset;
934 }
935
940 ~LanczosFunction1() noexcept override = default;
941
942 std::shared_ptr<Function1<ReturnT>> clone() const override {
944 }
945
946 ReturnT operator()(double x) const noexcept(IS_NOTHROW_INIT<ReturnT>) override {
947 double xArg1 = (x - this->_params[0]) * lsst::geom::PI;
948 double xArg2 = xArg1 * _invN;
949 if (std::fabs(xArg1) > 1.0e-5) {
950 return static_cast<ReturnT>(std::sin(xArg1) * std::sin(xArg2) / (xArg1 * xArg2));
951 } else {
952 return static_cast<ReturnT>(1);
953 }
954 }
955
959 unsigned int getOrder() const noexcept { return static_cast<unsigned int>(0.5 + (1.0 / _invN)); };
960
961 std::string toString(std::string const& prefix) const override {
963 os << "LanczosFunction1 [" << this->getOrder() << "]: ";
964 ;
965 os << Function1<ReturnT>::toString(prefix);
966 return os.str();
967 }
968
969private:
970 double _invN; // == 1/n
971
972protected:
973 /* Default constructor: intended only for serialization */
974 explicit LanczosFunction1() : Function1<ReturnT>(1), _invN(1.0) {}
975};
976
989template <typename ReturnT>
990class LanczosFunction2 : public Function2<ReturnT> {
991public:
995 explicit LanczosFunction2(unsigned int n,
996 double xOffset = 0.0,
997 double yOffset = 0.0)
998 : Function2<ReturnT>(2), _invN(1.0 / static_cast<double>(n)) {
999 this->_params[0] = xOffset;
1000 this->_params[1] = yOffset;
1001 }
1002
1007 ~LanczosFunction2() noexcept override = default;
1008
1009 std::shared_ptr<Function2<ReturnT>> clone() const override {
1011 new LanczosFunction2(this->getOrder(), this->_params[0], this->_params[1]));
1012 }
1013
1014 ReturnT operator()(double x, double y) const noexcept(IS_NOTHROW_INIT<ReturnT>) override {
1015 double xArg1 = (x - this->_params[0]) * lsst::geom::PI;
1016 double xArg2 = xArg1 * _invN;
1017 double xFunc = 1;
1018 if (std::fabs(xArg1) > 1.0e-5) {
1019 xFunc = std::sin(xArg1) * std::sin(xArg2) / (xArg1 * xArg2);
1020 }
1021 double yArg1 = (y - this->_params[1]) * lsst::geom::PI;
1022 double yArg2 = yArg1 * _invN;
1023 double yFunc = 1;
1024 if (std::fabs(yArg1) > 1.0e-5) {
1025 yFunc = std::sin(yArg1) * std::sin(yArg2) / (yArg1 * yArg2);
1026 }
1027 return static_cast<ReturnT>(xFunc * yFunc);
1028 }
1029
1033 unsigned int getOrder() const noexcept { return static_cast<unsigned int>(0.5 + (1.0 / _invN)); };
1034
1035 std::string toString(std::string const& prefix) const override {
1037 os << "LanczosFunction2 [" << this->getOrder() << "]: ";
1038 ;
1039 os << Function2<ReturnT>::toString(prefix);
1040 return os.str();
1041 }
1042
1043private:
1044 double _invN;
1045
1046protected:
1047 /* Default constructor: intended only for serialization */
1048 explicit LanczosFunction2() : Function2<ReturnT>(2), _invN(1.0) {}
1049};
1050} // namespace math
1051} // namespace afw
1052} // namespace lsst
1053
1054#endif // #ifndef LSST_AFW_MATH_FUNCTIONLIBRARY_H
#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
std::ostream * os
Definition Schema.cc:557
std::string prefix
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.
std::shared_ptr< Function1< ReturnT > > clone() const override
Return a pointer to a deep copy of this function.
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
Chebyshev1Function1 & operator=(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
Chebyshev1Function1 & operator=(Chebyshev1Function1 const &)=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.
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.
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
std::shared_ptr< Function2< ReturnT > > clone() const override
Return a pointer to a deep copy of this function.
~Chebyshev1Function2() noexcept override=default
lsst::geom::Box2D getXYRange() const noexcept
Get x,y range.
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 & operator=(Chebyshev1Function2 const &)=default
Chebyshev1Function2 & operator=(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 & operator=(DoubleGaussianFunction2 const &)=default
DoubleGaussianFunction2(double sigma1, double sigma2=0, double ampl2=0)
Construct a Gaussian function with specified x and y sigma.
DoubleGaussianFunction2(DoubleGaussianFunction2 &&)=default
DoubleGaussianFunction2(DoubleGaussianFunction2 const &)=default
std::string toString(std::string const &prefix) const override
Return a string representation of the function.
std::shared_ptr< Function2< ReturnT > > clone() const override
Return a pointer to a deep copy of this 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
std::shared_ptr< Function1< ReturnT > > clone() const override
Return a pointer to a deep copy of this function.
GaussianFunction1 & operator=(GaussianFunction1 &&)=default
GaussianFunction1(GaussianFunction1 &&)=default
std::string toString(std::string const &prefix) const override
Return a string representation of the function.
GaussianFunction1(GaussianFunction1 const &)=default
GaussianFunction1 & operator=(GaussianFunction1 const &)=default
GaussianFunction1(double sigma)
Construct a Gaussian function with specified sigma.
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 &&)=default
GaussianFunction2 & operator=(GaussianFunction2 const &)=default
std::shared_ptr< Function2< ReturnT > > clone() const override
Return a pointer to a deep copy of this function.
ReturnT operator()(double x, double y) const noexcept(IS_NOTHROW_INIT< ReturnT >) override
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(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.
IntegerDeltaFunction1 & operator=(IntegerDeltaFunction1 const &)=default
2-dimensional integer delta function.
IntegerDeltaFunction2 & operator=(IntegerDeltaFunction2 &&)=default
std::string toString(std::string const &prefix) const override
Return a string representation of the function.
IntegerDeltaFunction2(IntegerDeltaFunction2 const &)=default
ReturnT operator()(double x, double y) const noexcept(IS_NOTHROW_INIT< ReturnT >) override
IntegerDeltaFunction2 & operator=(IntegerDeltaFunction2 const &)=default
~IntegerDeltaFunction2() noexcept override=default
IntegerDeltaFunction2(double xo, double yo)
Construct an integer delta function with specified xo, yo.
std::shared_ptr< Function2< ReturnT > > clone() const override
Return a pointer to a deep copy of this function.
IntegerDeltaFunction2(IntegerDeltaFunction2 &&)=default
1-dimensional Lanczos function
LanczosFunction1(unsigned int n, double xOffset=0.0)
Construct a Lanczos function of specified order and x,y offset.
std::shared_ptr< Function1< ReturnT > > clone() const override
Return a pointer to a deep copy of this function.
LanczosFunction1 & operator=(LanczosFunction1 &&)=default
~LanczosFunction1() noexcept override=default
LanczosFunction1 & operator=(LanczosFunction1 const &)=default
LanczosFunction1(LanczosFunction1 &&)=default
LanczosFunction1(LanczosFunction1 const &)=default
ReturnT operator()(double x) const noexcept(IS_NOTHROW_INIT< ReturnT >) override
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
~LanczosFunction2() noexcept override=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
LanczosFunction2 & operator=(LanczosFunction2 const &)=default
unsigned int getOrder() const noexcept
Get the order of Lanczos function.
LanczosFunction2 & operator=(LanczosFunction2 &&)=default
std::string toString(std::string const &prefix) const override
Return a string representation of the function.
std::shared_ptr< Function2< ReturnT > > clone() const override
Return a pointer to a deep copy of this function.
1-dimensional polynomial function.
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.
bool isLinearCombination() const noexcept override
Is the function a linear combination of its parameters?
PolynomialFunction1 & operator=(PolynomialFunction1 &&)=default
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
std::shared_ptr< Function1< ReturnT > > clone() const override
Return a pointer to a deep copy of this function.
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 & operator=(PolynomialFunction2 &&)=default
std::shared_ptr< Function2< ReturnT > > clone() const override
Return a pointer to a deep copy of this function.
PolynomialFunction2(unsigned int order)
Construct a polynomial function of specified order.
std::string toString(std::string const &prefix) const override
Return a string representation of the function.
PolynomialFunction2(PolynomialFunction2 const &)=default
~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)
double constexpr PI
The ratio of a circle's circumference to diameter.
Definition Angle.h:40
STL namespace.
T sin(T... args)
T size(T... args)
table::Key< int > order