FunctionLibrary.h

Go to the documentation of this file.

// -*- LSST-C++ -*-

/*
 * LSST Data Management System
 * Copyright 2008, 2009, 2010 LSST Corporation.
 *
 * This product includes software developed by the
 * LSST Project (http://www.lsst.org/).
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the LSST License Statement and
 * the GNU General Public License along with this program.  If not,
 * see <http://www.lsstcorp.org/LegalNotices/>.
 */

#ifndef LSST_AFW_MATH_FUNCTIONLIBRARY_H
#define LSST_AFW_MATH_FUNCTIONLIBRARY_H

#include <algorithm>
#include <cmath>

#include "lsst/afw/geom.h"
#include "lsst/afw/math/Function.h"
#include "lsst/afw/geom/Angle.h"

namespace lsst {
namespace afw {
namespace math {

#ifndef SWIG
using boost::serialization::make_nvp;
#endif

    template<typename ReturnT>
    class IntegerDeltaFunction1: public Function1<ReturnT> {
    public:
        typedef typename Function1<ReturnT>::Ptr Function1Ptr;

        explicit IntegerDeltaFunction1(
            double xo)
        :
            Function1<ReturnT>(0),
            _xo(xo)
        {}

        virtual ~IntegerDeltaFunction1() {};

        virtual Function1Ptr clone() const {
            return Function1Ptr(new IntegerDeltaFunction1(_xo));
        }

        virtual ReturnT operator() (double x, double y) const {
            return static_cast<ReturnT>(x == _xo);
        }

        virtual std::string toString(std::string const& prefix="") const {
            std::ostringstream os;
            os << "IntegerDeltaFunction1 [" << _xo << "]: ";
            os << Function1<ReturnT>::toString(prefix);
            return os.str();
        };

    private:
        double _xo;

    protected:
        /* Default constructor: intended only for serialization */
        explicit IntegerDeltaFunction1() : Function1<ReturnT>(0), _xo(0.0) {}

    private:
        friend class boost::serialization::access;
        template <class Archive>
        void serialize(Archive& ar, unsigned int const version) {
            ar & make_nvp("fn1", boost::serialization::base_object<Function1<ReturnT> >(*this));
            ar & make_nvp("xo", this->_xo);
        }
    };

    template<typename ReturnT>
    class IntegerDeltaFunction2: public Function2<ReturnT> {
    public:
        typedef typename Function2<ReturnT>::Ptr Function2Ptr;

        explicit IntegerDeltaFunction2(
            double xo,
            double yo)
        :
            Function2<ReturnT>(0),
            _xo(xo),
            _yo(yo)
        {}

        virtual ~IntegerDeltaFunction2() {}

        virtual Function2Ptr clone() const {
            return Function2Ptr(new IntegerDeltaFunction2(_xo, _yo));
        }

        virtual ReturnT operator() (double x, double y) const {
            return static_cast<ReturnT>((x == _xo) && (y == _yo));
        }

        virtual std::string toString(std::string const& prefix) const {
            std::ostringstream os;
            os << "IntegerDeltaFunction2 [" << _xo << ", " << _yo << "]: ";
            os << Function2<ReturnT>::toString(prefix);
            return os.str();
        }

    private:
        double _xo;
        double _yo;

    protected:
        /* Default constructor: intended only for serialization */
        explicit IntegerDeltaFunction2() : Function2<ReturnT>(0), _xo(0.0), _yo(0.0) {}

    private:
        friend class boost::serialization::access;
        template <class Archive>
        void serialize(Archive& ar, unsigned int const version) {
            ar & make_nvp("fn2", boost::serialization::base_object<Function2<ReturnT> >(*this));
            ar & make_nvp("xo", this->_xo);
            ar & make_nvp("yo", this->_yo);
        }
    };

    template<typename ReturnT>
    class GaussianFunction1: public Function1<ReturnT> {
    public:
        typedef typename Function1<ReturnT>::Ptr Function1Ptr;

        explicit GaussianFunction1(
            double sigma)    
        :
            Function1<ReturnT>(1),
            _multFac(1.0 / std::sqrt(lsst::afw::geom::TWOPI))
        {
            this->_params[0] = sigma;
        }
        virtual ~GaussianFunction1() {}

        virtual Function1Ptr clone() const {
            return Function1Ptr(new GaussianFunction1(this->_params[0]));
        }

        virtual ReturnT operator() (double x) const {
            return static_cast<ReturnT> (
                (_multFac / this->_params[0]) *
                std::exp(- (x * x) / (2.0 * this->_params[0] * this->_params[0])));
        }

        virtual std::string toString(std::string const& prefix) const {
            std::ostringstream os;
            os << "GaussianFunction1 [" << _multFac << "]: ";
            os << Function1<ReturnT>::toString(prefix);
            return os.str();
        }

    private:
        const double _multFac; 

    protected:
        /* Default constructor: intended only for serialization */
        explicit GaussianFunction1() : Function1<ReturnT>(1), _multFac(1.0 / std::sqrt(lsst::afw::geom::TWOPI)) {}

    private:
        friend class boost::serialization::access;
        template <class Archive>
        void serialize(Archive& ar, unsigned int const version) {
            ar & make_nvp("fn1", boost::serialization::base_object<Function1<ReturnT> >(*this));
        }
    };

    template<typename ReturnT>
    class GaussianFunction2: public Function2<ReturnT> {
    public:
        typedef typename Function2<ReturnT>::Ptr Function2Ptr;

        explicit GaussianFunction2(
            double sigma1,      
            double sigma2,      
            double angle = 0.0) 
        :
            Function2<ReturnT>(3),
            _multFac(1.0 / (lsst::afw::geom::TWOPI))
        {
            this->_params[0] = sigma1;
            this->_params[1] = sigma2;
            this->_params[2] = angle;
            _updateCache();
        }

        virtual ~GaussianFunction2() {}

        virtual Function2Ptr clone() const {
            return Function2Ptr(new GaussianFunction2(this->_params[0], this->_params[1], this->_params[2]));
        }

        virtual ReturnT operator() (double x, double y) const {
            if (_angle != this->_params[2]) {
                _updateCache();
            }
            double pos1 = ( _cosAngle * x) + (_sinAngle * y);
            double pos2 = (-_sinAngle * x) + (_cosAngle * y);
            return static_cast<ReturnT> (
                (_multFac / (this->_params[0] * this->_params[1])) *
                std::exp(- ((pos1 * pos1) / (2.0 * this->_params[0] * this->_params[0]))
                         - ((pos2 * pos2) / (2.0 * this->_params[1] * this->_params[1]))));
        }

        virtual std::string toString(std::string const& prefix) const {
            std::ostringstream os;
            os << "GaussianFunction2: ";
            os << Function2<ReturnT>::toString(prefix);
            return os.str();
        }

        virtual bool isPersistable() const { return true; }

    protected:

        virtual std::string getPersistenceName() const;

        virtual void write(afw::table::io::OutputArchiveHandle & handle) const;

    private:
        void _updateCache() const {
            _angle = this->_params[2];
            _sinAngle = std::sin(_angle);
            _cosAngle = std::cos(_angle);
        }
        const double _multFac;  
        mutable double _angle;    
        mutable double _sinAngle; 
        mutable double _cosAngle; 

    protected:
        /* Default constructor: intended only for serialization */
        explicit GaussianFunction2() : Function2<ReturnT>(3), _multFac(1.0 / (lsst::afw::geom::TWOPI)), _angle(0.0),
            _sinAngle(0.0), _cosAngle(1.0) {}

    private:
        friend class boost::serialization::access;
        template <class Archive>
        void serialize(Archive& ar, unsigned int const version) {
            ar & make_nvp("fn2", boost::serialization::base_object<Function2<ReturnT> >(*this));
            ar & make_nvp("angle", this->_angle);
            ar & make_nvp("sinAngle", this->_sinAngle);
            ar & make_nvp("cosAngle", this->_cosAngle);
        }
    };

    template<typename ReturnT>
    class DoubleGaussianFunction2: public Function2<ReturnT> {
    public:
        typedef typename Function2<ReturnT>::Ptr Function2Ptr;

        explicit DoubleGaussianFunction2(
            double sigma1,      
            double sigma2 = 0,  
            double ampl2 = 0)   
        :
            Function2<ReturnT>(3),
            _multFac(1.0 / (lsst::afw::geom::TWOPI))
        {
            this->_params[0] = sigma1;
            this->_params[1] = sigma2;
            this->_params[2] = ampl2;
        }

        virtual ~DoubleGaussianFunction2() {}

        virtual Function2Ptr clone() const {
            return Function2Ptr(
                new DoubleGaussianFunction2(this->_params[0], this->_params[1], this->_params[2]));
        }

        virtual ReturnT operator() (double x, double y) const {
            double radSq = (x * x) + (y * y);
            double sigma1Sq = this->_params[0] * this->_params[0];
            double sigma2Sq = this->_params[1] * this->_params[1];
            double b = this->_params[2];
            return static_cast<ReturnT> (
                (_multFac / (sigma1Sq + (b * sigma2Sq))) *
                (std::exp(-radSq / (2.0 * sigma1Sq))
                + (b * std::exp(-radSq / (2.0 * sigma2Sq)))));
        }

        virtual std::string toString(std::string const& prefix) const {
            std::ostringstream os;
            os << "DoubleGaussianFunction2 [" << _multFac << "]: ";
            os << Function2<ReturnT>::toString(prefix);
            return os.str();
        }

        virtual bool isPersistable() const { return true; }

    protected:

        virtual std::string getPersistenceName() const;

        virtual void write(afw::table::io::OutputArchiveHandle & handle) const;

    private:
        const double _multFac; 

    protected:
        /* Default constructor: intended only for serialization */
        explicit DoubleGaussianFunction2() : Function2<ReturnT>(3), _multFac(1.0 / (lsst::afw::geom::TWOPI)) {}

    private:
        friend class boost::serialization::access;
        template <class Archive>
        void serialize(Archive& ar, unsigned int const version) {
            ar & make_nvp("fn2", boost::serialization::base_object<Function2<ReturnT> >(*this));
        }
    };

    template<typename ReturnT>
    class PolynomialFunction1: public Function1<ReturnT> {
    public:
        typedef typename Function1<ReturnT>::Ptr Function1Ptr;

        explicit PolynomialFunction1(
            unsigned int order)     
        :
            Function1<ReturnT>(order+1) {
        }

        explicit PolynomialFunction1(
            std::vector<double> params)  
        :
            Function1<ReturnT>(params)
        {
            if (params.size() < 1) {
                throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterError,
                                  "PolynomialFunction1 called with empty vector");
            }
        }

        virtual ~PolynomialFunction1() {}

        virtual Function1Ptr clone() const {
            return Function1Ptr(new PolynomialFunction1(this->_params));
        }

        virtual bool isLinearCombination() const { return true; };

        virtual ReturnT operator() (double x) const {
            int const order = static_cast<int>(this->_params.size()) - 1;
            double retVal = this->_params[order];
            for (int ii = order-1; ii >= 0; --ii) {
                retVal = (retVal * x) + this->_params[ii];
            }
            return static_cast<ReturnT>(retVal);
        }

        unsigned int getOrder() const { return this->getNParameters() - 1; };

        virtual std::string toString(std::string const& prefix) const {
            std::ostringstream os;
            os << "PolynomialFunction1 []: ";
            os << Function1<ReturnT>::toString(prefix);
            return os.str();
        }

    protected:
        /* Default constructor: intended only for serialization */
        explicit PolynomialFunction1() : Function1<ReturnT>(1) {}

    private:
        friend class boost::serialization::access;
        template <class Archive>
        void serialize(Archive& ar, unsigned int const version) {
            ar & make_nvp("fn1", boost::serialization::base_object<Function1<ReturnT> >(*this));
        }
    };

    template<typename ReturnT>
    class PolynomialFunction2: public BasePolynomialFunction2<ReturnT> {
    public:
        typedef typename Function2<ReturnT>::Ptr Function2Ptr;

        explicit PolynomialFunction2(
            unsigned int order) 
        :
            BasePolynomialFunction2<ReturnT>(order),
            _oldY(0),
            _xCoeffs(this->_order + 1)
        {}

        explicit PolynomialFunction2(
            std::vector<double> params)  
        :
            BasePolynomialFunction2<ReturnT>(params),
            _oldY(0),
            _xCoeffs(this->_order + 1)
        {}

        virtual ~PolynomialFunction2() {}

        virtual Function2Ptr clone() const {
            return Function2Ptr(new PolynomialFunction2(this->_params));
        }

        virtual ReturnT operator() (double x, double y) const {
            /* Solve as follows:
            - f(x,y) = Cx0 + Cx1 x + Cx2 x^2 + Cx3 x^3 + ...
            where:
              Cx0 = P0 + P2 y + P5 y^2 + P9 y^3 + ...
              Cx1 = P1 + P4 y + P8 y2 + ...
              Cx2 = P3 + P7 y + ...
              Cx3 = P6 + ...
              ...

            Compute Cx0, Cx1...Cxn by solving 1-d polynomials in y in the usual way.
            These values are cached and only recomputed for new values of Y or if the parameters change.

            Then compute f(x,y) by solving the 1-d polynomial in x in the usual way.
            */
            const int maxXCoeffInd = this->_order;

            if ((y != _oldY) || !this->_isCacheValid) {
                // update _xCoeffs cache
                // note: paramInd is decremented in both of the following loops
                int paramInd = static_cast<int>(this->_params.size()) - 1;

                // initialize _xCoeffs to coeffs for pure y^n; e.g. for 3rd order:
                // _xCoeffs[0] = _params[9], _xCoeffs[1] = _params[8], ... _xCoeffs[3] = _params[6]
                for (int xCoeffInd = 0; xCoeffInd <= maxXCoeffInd; ++xCoeffInd, --paramInd) {
                    _xCoeffs[xCoeffInd] = this->_params[paramInd];
                }

                // finish computing _xCoeffs
                for (int xCoeffInd = 0, endXCoeffInd = maxXCoeffInd; paramInd >= 0; --paramInd) {
                    _xCoeffs[xCoeffInd] = (_xCoeffs[xCoeffInd] * y) + this->_params[paramInd];
                    ++xCoeffInd;
                    if (xCoeffInd >= endXCoeffInd) {
                        xCoeffInd = 0;
                        --endXCoeffInd;
                    }
                }

                _oldY = y;
                this->_isCacheValid = true;
            }

            // use _xCoeffs to compute result
            double retVal = _xCoeffs[maxXCoeffInd];
            for (int xCoeffInd = maxXCoeffInd - 1; xCoeffInd >= 0; --xCoeffInd) {
                retVal = (retVal * x) + _xCoeffs[xCoeffInd];
            }
            return static_cast<ReturnT>(retVal);
        }

        virtual std::vector<double> getDFuncDParameters(double x, double y) const;

        virtual std::string toString(std::string const& prefix) const {
            std::ostringstream os;
            os << "PolynomialFunction2 [" << this->_order << "]: ";
            os << Function2<ReturnT>::toString(prefix);
            return os.str();
        }

        virtual bool isPersistable() const { return true; }

    protected:

        virtual std::string getPersistenceName() const;

        virtual void write(afw::table::io::OutputArchiveHandle & handle) const;

    private:
        mutable double _oldY;         
        mutable std::vector<double> _xCoeffs; 

    protected:
        /* Default constructor: intended only for serialization */
        explicit PolynomialFunction2() : BasePolynomialFunction2<ReturnT>(), _oldY(0), _xCoeffs(0)  {}

    private:
        friend class boost::serialization::access;
        template <class Archive>
        void serialize(Archive& ar, unsigned int const version) {
            ar & make_nvp("fn2", boost::serialization::base_object<BasePolynomialFunction2<ReturnT> >(*this));
            ar & make_nvp("yCoeffs", this->_xCoeffs); // sets size of _xCoeffs; name is historical
        }
    };

    template<typename ReturnT>
    class Chebyshev1Function1: public Function1<ReturnT> {
    public:
        typedef typename Function1<ReturnT>::Ptr Function1Ptr;

        explicit Chebyshev1Function1(
            unsigned int order, 
            double minX = -1,   
            double maxX = 1)    
        :
            Function1<ReturnT>(order + 1)
        {
            _initialize(minX, maxX);
        }

        explicit Chebyshev1Function1(
            std::vector<double> params, 
            double minX = -1,   
            double maxX = 1)    
        :
            Function1<ReturnT>(params)
        {
            if (params.size() < 1) {
                throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterError,
                                  "Chebyshev1Function1 called with empty vector");
            }
            _initialize(minX, maxX);
        }

        virtual ~Chebyshev1Function1() {}

        virtual Function1Ptr clone() const {
            return Function1Ptr(new Chebyshev1Function1(this->_params, _minX, _maxX));
        }

        double getMinX() const { return _minX; };

        double getMaxX() const { return _maxX; };

        unsigned int getOrder() const { return this->getNParameters() - 1; };

        virtual bool isLinearCombination() const { return true; };

        virtual ReturnT operator() (double x) const {
            double xPrime = (x + _offset) * _scale;
            
            // Clenshaw function for solving the Chebyshev polynomial
            // Non-recursive version from Kresimir Cosic
            int const order = _order;
            if (order == 0) {
                return this->_params[0];
            } else if (order == 1) {
                return this->_params[0] + (this->_params[1] * xPrime);
            }
            double cshPrev = this->_params[order];
            double csh = (2 * xPrime * this->_params[order]) + this->_params[order-1];
            for (int i = order - 2; i > 0; --i) {
                double cshNext = (2 * xPrime * csh) + this->_params[i] - cshPrev;
                cshPrev = csh;
                csh = cshNext;
            }
            return (xPrime * csh) + this->_params[0] - cshPrev;
        }

        virtual std::string toString(std::string const& prefix) const {
            std::ostringstream os;
            os << "Chebyshev1Function1 [" << _minX << ", " << _maxX << "]: ";
            os << Function1<ReturnT>::toString(prefix);
            return os.str();
        }

    private:
        double _minX;    
        double _maxX;    
        double _scale;   
        double _offset;  
        unsigned int _order;   

        void _initialize(double minX, double maxX) {
            _minX = minX;
            _maxX = maxX;
            _scale = 2.0 / (_maxX - _minX);
            _offset = -(_minX + _maxX) * 0.5;
            _order = this->getNParameters() - 1;
        }

    protected:
        /* Default constructor: intended only for serialization */
        explicit Chebyshev1Function1() : Function1<ReturnT>(1),
            _minX(0.0), _maxX(0.0), _scale(1.0), _offset(0.0), _order(0) {}

    private:
        friend class boost::serialization::access;
        template <class Archive>
        void serialize(Archive& ar, unsigned int const) {
            ar & make_nvp("fn1", boost::serialization::base_object<Function1<ReturnT> >(*this));
            ar & make_nvp("minX", this->_minX);
            ar & make_nvp("minX", this->_minX);
            ar & make_nvp("maxX", this->_maxX);
            ar & make_nvp("scale", this->_scale);
            ar & make_nvp("offset", this->_offset);
            ar & make_nvp("maxInd", this->_order);
        }
    };

    template<typename ReturnT>
    class Chebyshev1Function2: public BasePolynomialFunction2<ReturnT> {
    public:
        typedef typename Function2<ReturnT>::Ptr Function2Ptr;

        explicit Chebyshev1Function2(
            unsigned int order, 
            lsst::afw::geom::Box2D const &xyRange =
                lsst::afw::geom::Box2D(lsst::afw::geom::Point2D(-1.0, -1.0),
                                       lsst::afw::geom::Point2D( 1.0,  1.0)))   
        :
            BasePolynomialFunction2<ReturnT>(order),
            _oldYPrime(0),
            _yCheby(this->_order + 1),
            _xCoeffs(this->_order + 1)
        {
            _initialize(xyRange);
        }

        explicit Chebyshev1Function2(
            std::vector<double> params, 
            lsst::afw::geom::Box2D const &xyRange =
                lsst::afw::geom::Box2D(lsst::afw::geom::Point2D(-1.0, -1.0),
                                       lsst::afw::geom::Point2D( 1.0,  1.0)))   
        :
            BasePolynomialFunction2<ReturnT>(params),
            _oldYPrime(0),
            _yCheby(this->_order + 1),
            _xCoeffs(this->_order + 1)
        {
            _initialize(xyRange);
        }

        virtual ~Chebyshev1Function2() {}

        virtual Function2Ptr clone() const {
            return Function2Ptr(new Chebyshev1Function2(this->_params, this->getXYRange()));
        }

        lsst::afw::geom::Box2D getXYRange() const {
            return lsst::afw::geom::Box2D(lsst::afw::geom::Point2D(_minX, _minY),
                                          lsst::afw::geom::Point2D(_maxX, _maxY));
        };

        virtual Chebyshev1Function2 truncate(
                int truncOrder 
        ) const {
            if (truncOrder > this->_order) {
                std::ostringstream os;
                os << "truncated order=" << truncOrder << " must be <= original order=" << this->_order;
                throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterError, os.str());
            }
            int truncNParams = this->nParametersFromOrder(truncOrder);
            std::vector<double> truncParams(this->_params.begin(), this->_params.begin() + truncNParams);
            return Chebyshev1Function2(truncParams, this->getXYRange());
        }

        virtual ReturnT operator() (double x, double y) const {
            /* Solve as follows:
            - f(x,y) = Cy0 T0(y') + Cy1 T1(y') + Cy2 T2(y') + Cy3 T3(y') + ...
            where:
              Cy0 = P0 T0(x') + P1 T1(x') + P3 T2(x') + P6 T3(x') + ...
              Cy1 = P2 T0(x') + P4 T1(x') + P7 T2(x') + ...
              Cy2 = P5 T0(x') + P8 T1(x') + ...
              Cy3 = P9 T0(x') + ...
              ...

            First compute Tn(x') for each n
            Then use that to compute Cy0, Cy1, ...Cyn
            Then solve the y Chebyshev polynomial using the Clenshaw algorithm
            */
            double const xPrime = (x + _offsetX) * _scaleX;
            double const yPrime = (y + _offsetY) * _scaleY;

            const int nParams = static_cast<int>(this->_params.size());
            const int order = this->_order;

            if (order == 0) {
                return this->_params[0];  // No caching required
            }

            if ((yPrime != _oldYPrime) || !this->_isCacheValid) {
                // update cached _yCheby and _xCoeffs
                _yCheby[0] = 1.0;
                _yCheby[1] = yPrime;
                for (int chebyInd = 2; chebyInd <= order; chebyInd++) {
                    _yCheby[chebyInd] = (2 * yPrime * _yCheby[chebyInd-1]) - _yCheby[chebyInd-2];
                }

                for (int coeffInd=0; coeffInd <= order; coeffInd++) {
                    _xCoeffs[coeffInd] = 0;
                }
                for (int coeffInd = 0, endCoeffInd = 0, paramInd = 0; paramInd < nParams; paramInd++) {
                    _xCoeffs[coeffInd] += this->_params[paramInd] * _yCheby[endCoeffInd];
                    --coeffInd;
                    ++endCoeffInd;
                    if (coeffInd < 0) {
                        coeffInd = endCoeffInd;
                        endCoeffInd = 0;
                    }
                }

                _oldYPrime = yPrime;
                this->_isCacheValid = true;
            }

            // Clenshaw function for solving the Chebyshev polynomial
            // Non-recursive version from Kresimir Cosic
            if (order == 1) {
                return _xCoeffs[0] + (_xCoeffs[1] * xPrime);
            }
            double cshPrev = _xCoeffs[order];
            double csh = (2 * xPrime * _xCoeffs[order]) + _xCoeffs[order-1];
            for (int i = order - 2; i > 0; --i) {
                double cshNext = (2 * xPrime * csh) + _xCoeffs[i] - cshPrev;
                cshPrev = csh;
                csh = cshNext;
            }
            return (xPrime * csh) + _xCoeffs[0] - cshPrev;
        }

        virtual std::string toString(std::string const& prefix) const {
            std::ostringstream os;
            os << "Chebyshev1Function2 [";
            os << this->_order << ", " << this->getXYRange() << "]";
            os << Function2<ReturnT>::toString(prefix);
            return os.str();
        }

        virtual bool isPersistable() const { return true; }

    protected:

        virtual std::string getPersistenceName() const;

        virtual void write(afw::table::io::OutputArchiveHandle & handle) const;

    private:
        mutable double _oldYPrime;
        mutable std::vector<double> _yCheby;    
        mutable std::vector<double> _xCoeffs;   
        double _minX;    
        double _minY;    
        double _maxX;    
        double _maxY;    
        double _scaleX;   
        double _scaleY;   
        double _offsetX;  
        double _offsetY;  

        void _initialize(lsst::afw::geom::Box2D const &xyRange) {
            _minX = xyRange.getMinX();
            _minY = xyRange.getMinY();
            _maxX = xyRange.getMaxX();
            _maxY = xyRange.getMaxY();
            _scaleX = 2.0 / (_maxX - _minX);
            _scaleY = 2.0 / (_maxY - _minY);
            _offsetX = -(_minX + _maxX) * 0.5;
            _offsetY = -(_minY + _maxY) * 0.5;
        }

    protected:
        /* Default constructor: intended only for serialization */
        explicit Chebyshev1Function2() : BasePolynomialFunction2<ReturnT>(),
            _oldYPrime(0),
            _yCheby(0),
            _xCoeffs(0),
            _minX(0.0), _minY(0.0),
            _maxX(0.0), _maxY(0.0),
            _scaleX(1.0), _scaleY(1.0),
            _offsetX(0.0), _offsetY(0.0) {}

    private:
        friend class boost::serialization::access;
        template <class Archive>
        void serialize(Archive& ar, unsigned int const version) {
            ar & make_nvp("fn2", boost::serialization::base_object<BasePolynomialFunction2<ReturnT> >(*this));
            ar & make_nvp("minX", this->_minX);
            ar & make_nvp("minY", this->_minY);
            ar & make_nvp("maxX", this->_maxX);
            ar & make_nvp("maxY", this->_maxY);
            ar & make_nvp("scaleX", this->_scaleX);
            ar & make_nvp("scaleY", this->_scaleY);
            ar & make_nvp("offsetX", this->_offsetX);
            ar & make_nvp("offsetY", this->_offsetY);
            ar & make_nvp("xCheby", this->_yCheby); // sets size of _yCheby; name is historical
            _xCoeffs.resize(_yCheby.size());
        }
    };

    template<typename ReturnT>
    class LanczosFunction1: public Function1<ReturnT> {
    public:
        typedef typename Function1<ReturnT>::Ptr Function1Ptr;

        explicit LanczosFunction1(
            unsigned int n,         
            double xOffset = 0.0)    
        :
            Function1<ReturnT>(1),
            _invN(1.0/static_cast<double>(n))
        {
            this->_params[0] = xOffset;
        }

        virtual ~LanczosFunction1() {}

        virtual Function1Ptr clone() const {
            return Function1Ptr(new LanczosFunction1(this->getOrder(), this->_params[0]));
        }

        virtual ReturnT operator() (double x) const {
            double xArg1 = (x - this->_params[0]) * lsst::afw::geom::PI;
            double xArg2 = xArg1 * _invN;
            if (std::fabs(xArg1) > 1.0e-5) {
                return static_cast<ReturnT>(std::sin(xArg1) * std::sin(xArg2) / (xArg1 * xArg2));
            } else {
                return static_cast<ReturnT>(1);
            }
        }

        unsigned int getOrder() const {
            return static_cast<unsigned int>(0.5 + (1.0 / _invN));
        };

        virtual std::string toString(std::string const& prefix) const {
            std::ostringstream os;
            os << "LanczosFunction1 [" << this->getOrder() << "]: ";;
            os << Function1<ReturnT>::toString(prefix);
            return os.str();
        }

    private:
        double _invN;                   // == 1/n

    protected:
        /* Default constructor: intended only for serialization */
        explicit LanczosFunction1() : Function1<ReturnT>(1), _invN(1.0) {}

    private:
        friend class boost::serialization::access;
        template <class Archive>
        void serialize(Archive& ar, unsigned int const version) {
            ar & make_nvp("fn1", boost::serialization::base_object<Function1<ReturnT> >(*this));
            ar & make_nvp("invN", this->_invN);
        }
    };

    template<typename ReturnT>
    class LanczosFunction2: public Function2<ReturnT> {
    public:
        typedef typename Function2<ReturnT>::Ptr Function2Ptr;

        explicit LanczosFunction2(
            unsigned int n,         
            double xOffset = 0.0,    
            double yOffset = 0.0)    
        :
            Function2<ReturnT>(2),
            _invN(1.0 / static_cast<double>(n))
        {
            this->_params[0] = xOffset;
            this->_params[1] = yOffset;
        }

        virtual ~LanczosFunction2() {}

        virtual Function2Ptr clone() const {
            return Function2Ptr(new LanczosFunction2(this->getOrder(), this->_params[0], this->_params[1]));
        }

        virtual ReturnT operator() (double x, double y) const {
            double xArg1 = (x - this->_params[0]) * lsst::afw::geom::PI;
            double xArg2 = xArg1 * _invN;
            double xFunc = 1;
            if (std::fabs(xArg1) > 1.0e-5) {
                xFunc = std::sin(xArg1) * std::sin(xArg2) / (xArg1 * xArg2);
            }
            double yArg1 = (y - this->_params[1]) * lsst::afw::geom::PI;
            double yArg2 = yArg1 * _invN;
            double yFunc = 1;
            if (std::fabs(yArg1) > 1.0e-5) {
                yFunc = std::sin(yArg1) * std::sin(yArg2) / (yArg1 * yArg2);
            }
            return static_cast<ReturnT>(xFunc * yFunc);
        }

        unsigned int getOrder() const {
            return static_cast<unsigned int>(0.5 + (1.0 / _invN));
        };

        virtual std::string toString(std::string const& prefix) const {
            std::ostringstream os;
            os << "LanczosFunction2 [" << this->getOrder() << "]: ";;
            os << Function2<ReturnT>::toString(prefix);
            return os.str();
        }

    private:
        double _invN;   

    protected:
        /* Default constructor: intended only for serialization */
        explicit LanczosFunction2() : Function2<ReturnT>(2), _invN(1.0) {}

    private:
        friend class boost::serialization::access;
        template <class Archive>
        void serialize(Archive& ar, unsigned int const version) {
            ar & make_nvp("fn2", boost::serialization::base_object<Function2<ReturnT> >(*this));
            ar & make_nvp("invN", this->_invN);
        }
    };

}}}   // lsst::afw::math

#endif // #ifndef LSST_AFW_MATH_FUNCTIONLIBRARY_H