lsst::afw::math::detail::TautSpline Class Reference

#include <Spline.h>

Inheritance diagram for lsst::afw::math::detail::TautSpline:

Public Types
enum	Symmetry { Unknown, Odd, Even }

Public Member Functions
	TautSpline (std::vector< double > const &x, std::vector< double > const &y, double const gamma=0, Symmetry type=Unknown)
Public Member Functions inherited from lsst::afw::math::detail::Spline
virtual	~Spline ()
void	interpolate (std::vector< double > const &x, std::vector< double > &y) const
void	derivative (std::vector< double > const &x, std::vector< double > &dydx) const
std::vector< double >	roots (double const value, double const x0, double const x1) const

Private Member Functions
void	calculateTautSpline (std::vector< double > const &x, std::vector< double > const &y, double const gamma0)
void	calculateTautSplineEvenOdd (std::vector< double > const &x, std::vector< double > const &y, double const gamma0, bool even)

Additional Inherited Members
Protected Member Functions inherited from lsst::afw::math::detail::Spline
	Spline ()
void	_allocateSpline (int const nknot)
Protected Attributes inherited from lsst::afw::math::detail::Spline
std::vector< double >	_knots
std::vector< std::vector < double > >	_coeffs

Detailed Description

Definition at line 36 of file Spline.h.

Member Enumeration Documentation

enum lsst::afw::math::detail::TautSpline::Symmetry

Enumerator
Unknown
Odd
Even

Definition at line 38 of file Spline.h.

{ Unknown, Odd, Even };

Constructor & Destructor Documentation

lsst::afw::math::detail::TautSpline::TautSpline	(	std::vector< double > const &	x,
		std::vector< double > const &	y,
		double const	gamma0 = 0,
		Symmetry	type = Unknown
	)

Adapted from A Practical Guide to Splines by C. de Boor (N.Y. : Springer-Verlag, 1978) (His routine tautsp converted to C by Robert Lupton and then to C++ by an older and grayer Robert Lupton)

Construct cubic spline interpolant to given data.

If gamma .gt. 0., additional knots are introduced where needed to make the interpolant more flexible locally. this avoids extraneous inflection points typical of cubic spline interpolation at knots to rapidly changing data. Values for gamma are:

     = 0., no additional knots
     in (0.,3.), under certain conditions on the given data at
           points i-1, i, i+1, and i+2, a knot is added in the
           i-th interval, i=1,...,ntau-3. see description of method
           below. the interpolant gets rounded with increasing
           gamma. a value of  2.5  for gamma is typical.
     in (3.,6.), same , except that knots might also be added in
           intervals in which an inflection point would be permitted.
           a value of  5.5  for gamma is typical.

Returns: a SPLINE giving knots and coefficients, if all was OK NULL otherwise (and pushes message on error stack)

---— m e t h o d ---— on the i-th interval, (tau[i], tau[i+1]), the interpolant is of the form (*) f(u(x)) = a + b*u + c*h(u,z) + d*h(1-u,1-z) , with u = u(x) = (x - tau[i])/dtau[i]. Here, z = z(i) = addg(i+1)/(addg(i) + addg(i+1)) (= .5, in case the denominator vanishes). with ddg(j) = dg(j+1) - dg(j), addg(j) = abs(ddg(j)), dg(j) = divdif(j) = (gtau[j+1] - gtau[j])/dtau[j] and h(u,z) = alpha*u**3 + (1 - alpha)*(max(((u-zeta)/(1-zeta)),0)**3 with alpha(z) = (1-gamma/3)/zeta zeta(z) = 1 - gamma*min((1 - z), 1/3) thus, for 1/3 .le. z .le. 2/3, f is just a cubic polynomial on the interval i. otherwise, it has one additional knot, at tau[i] + zeta*dtau[i] . as z approaches 1, h(.,z) has an increasingly sharp bend near 1, thus allowing f to turn rapidly near the additional knot. in terms of f(j) = gtau[j] and fsecnd[j] = 2*derivative of f at tau[j], the coefficients for (*) are given as a = f(i) - d b = (f(i+1) - f(i)) - (c - d) c = fsecnd[i+1]*dtau[i]**2/hsecnd(1,z) d = fsecnd[i]*dtau[i]**2/hsecnd(1,1-z) hence can be computed once fsecnd[i],i=0,...,ntau-1, is fixed.

f is automatically continuous and has a continuous second derivat- ive (except when z = 0 or 1 for some i). we determine fscnd(.) from the requirement that also the first derivative of f be continuous. in addition, we require that the third derivative be continuous across tau[1] and across tau[ntau-2] . this leads to a strictly diagonally dominant tridiagonal linear system for the fsecnd[i] which we solve by gauss elimination without pivoting.

There must be at least 4 interpolation points for us to fit a taut cubic spline, but if you provide fewer we'll fit a quadratic or linear polynomial (but you must provide at least 2)

See also phSplineFindTautEven and phSplineFindTautOdd for splining functions with known symmetries

Parameters

x	points where function's specified
y	values of function at tau[]
gamma0	control extra knots
type	specify the desired symmetry (e.g. Even)

Definition at line 178 of file Spline.cc.

                {
                    if(x.size() != y.size()) {
                        throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterError,
                                          (boost::format("TautSpline: x and y must have the same size; saw %d %d\n")
                                           % x.size() % y.size()).str());
                    }

                    int const ntau = x.size();          /* size of tau and gtau, must be >= 2*/
                    if(ntau < 2) {
                        throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterError,
                                          (boost::format("TautSpline: ntau = %d, should be >= 2\n")
                                           % ntau).str());
                    }

                    switch (type) {
                      case Unknown:
                        calculateTautSpline(x, y, gamma0);
                        break;
                      case Even:
                      case Odd:
                        calculateTautSplineEvenOdd(x, y, gamma0, type == Even);
                        break;
                    }
                }

Member Function Documentation

void lsst::afw::math::detail::TautSpline::calculateTautSpline ( std::vector< double > const &  x, std::vector< double > const &  y, double const  gamma0  )

private

Here's the worker routine for the TautSpline ctor

Parameters

x	points where function's specified
y	values of function at tau[]
gamma0	control extra knots

Definition at line 213 of file Spline.cc.

                {
                    const double *tau = &x[0];
                    const double *gtau = &y[0];
                    int const ntau = x.size();          // size of tau and gtau, must be >= 2

                    if(ntau < 4) {                      // use a single quadratic
                        int const nknot = ntau;

                        _allocateSpline(nknot);
        
                        _knots[0] = tau[0];
                        for (int i = 1; i < nknot;i++) {
                            _knots[i] = tau[i];
                            if(tau[i - 1] >= tau[i]) {
                                throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterError,
                                                  (boost::format("point %d and the next, %f %f, are out of order")
                                                   % (i - 1)  % tau[i - 1] % tau[i]).str());
                            }
                        }
        
                        if(ntau == 2) {
                            _coeffs[0][0] = gtau[0];
                            _coeffs[1][0] = (gtau[1] - gtau[0])/(tau[1] - tau[0]);
                            _coeffs[2][0] = _coeffs[3][0] = 0;

                            _coeffs[0][1] = gtau[1];
                            _coeffs[1][1] = (gtau[1] - gtau[0])/(tau[1] - tau[0]);
                            _coeffs[2][1] = _coeffs[3][1] = 0;
                        } else {                                /* must be 3 */
                            double tmp = (tau[2]-tau[0])*(tau[2]-tau[1])*(tau[1]-tau[0]);
                            _coeffs[0][0] = gtau[0];
                            _coeffs[1][0] = ((gtau[1] - gtau[0])*pow(tau[2] - tau[0],2) -
                                             (gtau[2] - gtau[0])*pow(tau[1] - tau[0],2))/tmp;
                            _coeffs[2][0] = -2*((gtau[1] - gtau[0])*(tau[2] - tau[0]) -
                                                (gtau[2] - gtau[0])*(tau[1] - tau[0]))/tmp;
                            _coeffs[3][0] = 0;

                            _coeffs[0][1] = gtau[1];
                            _coeffs[1][1] = _coeffs[1][0] + (tau[1] - tau[0])*_coeffs[2][0];
                            _coeffs[2][1] = _coeffs[2][0];
                            _coeffs[3][1] = 0;

                            _coeffs[0][2] = gtau[2];
                            _coeffs[1][2] = _coeffs[1][0] + (tau[2] - tau[0])*_coeffs[2][0];
                            _coeffs[2][2] = _coeffs[2][0];
                            _coeffs[3][2] = 0;
                        }
       
                        return;
                    }
/*
 * Allocate scratch space
 *     s[0][...] = dtau = tau(.+1) - tau
 *     s[1][...] = diag = diagonal in linear system
 *     s[2][...] = u = upper diagonal in linear system
 *     s[3][...] = r = right side for linear system (initially)
 *               = fsecnd = solution of linear system , namely the second
 *                          derivatives of interpolant at  tau
 *     s[4][...] = z = indicator of additional knots
 *     s[5][...] = 1/hsecnd(1,x) with x = z or = 1-z. see below.
 */
                    std::vector<std::vector<double> > s(6); // scratch space

                    for (int i = 0; i != 6; i++) {
                        s[i].resize(ntau);
                    }
/*
 * Construct delta tau and first and second (divided) differences of data 
 */

                    for (int i = 0; i < ntau - 1; i++) {
                        s[0][i] = tau[i + 1] - tau[i];
                        if(s[0][i] <= 0.) {
                            throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterError,
                                              (boost::format("point %d and the next, %f %f, are out of order")
                                               % i % tau[i] % tau[i+1]).str());
                        }
                        s[3][i + 1] = (gtau[i + 1] - gtau[i])/s[0][i];
                    }
                    for (int i = 1; i < ntau - 1; ++i) {
                        s[3][i] = s[3][i + 1] - s[3][i];
                    }
/*
 * Construct system of equations for second derivatives at tau. At each
 * interior data point, there is one continuity equation, at the first
 * and the last interior data point there is an additional one for a
 * total of ntau equations in ntau unknowns.
 */
                    s[1][1] = s[0][0]/3;

                    int method;
                    double gamma = gamma0;               // control the smoothing
                    if(gamma <= 0) {
                        method = 1;
                    } else if(gamma > 3) {
                        gamma -= 3;
                        method = 3;
                    } else {
                        method = 2;
                    }
                    double const onemg3 = 1 - gamma/3;

                    int nknot = ntau;                   // count number of knots
/*
 * Some compilers don't realise that the flow of control always initialises
 * these variables; so initialise them to placate e.g. gcc
 */
                    double zeta = 0.0;
                    double alpha = 0.0;
                    double ratio = 0.0;

                    //double c, d;
                    //double z, denom, factr2;
                    //double onemzt, zt2, del;


                    double entry3 = 0.0;
                    double factor = 0.0;
                    double onemzt = 0;
                    double zt2 = 0;
                    double z_half = 0;
    
                    for (int i = 1; i < ntau - 2; ++i) {
                        /*
                         * construct z[i] and zeta[i]
                         */
                        double z = .5;
                        if((method == 2 && s[3][i]*s[3][i + 1] >= 0) || method == 3) {
                            double const temp = fabs(s[3][i + 1]);
                            double const denom = fabs(s[3][i]) + temp;
                            if(denom != 0) {
                                z = temp/denom;
                                if (fabs(z - 0.5) <= 1.0/6.0) {
                                    z = 0.5;
                                }
                            }
                        }
       
                        s[4][i] = z;
/*
  set up part of the i-th equation which depends on the i-th interval
*/
                        z_half = z - 0.5;
                        if(z_half < 0) {
                            zeta = gamma*z;
                            onemzt = 1 - zeta;
                            zt2 = zeta*zeta;
                            double temp = onemg3/onemzt;
                            alpha = (temp < 1 ? temp : 1);
                            factor = zeta/(alpha*(zt2 - 1) + 1);
                            s[5][i] = zeta*factor/6;
                            s[1][i] += s[0][i]*((1 - alpha*onemzt)*factor/2 - s[5][i]);
/*
 * if z = 0 and the previous z = 1, then d[i] = 0. Since then
 * also u[i-1] = l[i+1] = 0, its value does not matter. Reset
 * d[i] = 1 to insure nonzero pivot in elimination.
 */
                            if (s[1][i] <= 0.) {
                                s[1][i] = 1;
                            }
                            s[2][i] = s[0][i]/6;
           
                            if(z != 0) {                        /* we'll get a new knot */
                                nknot++;
                            }
                        } else if(z_half == 0) {
                            s[1][i] += s[0][i]/3;
                            s[2][i] = s[0][i]/6;
                        } else {
                            onemzt = gamma*(1 - z);
                            zeta = 1 - onemzt;
                            double const temp = onemg3/zeta;
                            alpha = (temp < 1 ? temp : 1);
                            factor = onemzt/(1 - alpha*zeta*(onemzt + 1));
                            s[5][i] = onemzt*factor/6;
                            s[1][i] += s[0][i]/3;
                            s[2][i] = s[5][i]*s[0][i];
           
                            if(onemzt != 0) {                   /* we'll get a new knot */
                                nknot++;
                            }
                        }
                        if (i == 1) {
                            s[4][0] = 0.5;
/*
 * the first two equations enforce continuity of the first and of
 * the third derivative across tau[1].
 */
                            s[1][0] = s[0][0]/6;
                            s[2][0] = s[1][1];
                            entry3 = s[2][1];
                            if(z_half < 0) {
                                const double factr2 = zeta*(alpha*(zt2 - 1.) + 1.)/(alpha*(zeta*zt2 - 1.) + 1.);
                                ratio = factr2*s[0][1]/s[1][0];
                                s[1][1] = factr2*s[0][1] + s[0][0];
                                s[2][1] = -factr2*s[0][0];
                            } else if (z_half == 0) {
                                ratio = s[0][1]/s[1][0];
                                s[1][1] = s[0][1] + s[0][0];
                                s[2][1] = -s[0][0];
                            } else {
                                ratio = s[0][1]/s[1][0];
                                s[1][1] = s[0][1] + s[0][0];
                                s[2][1] = -s[0][0]*6*alpha*s[5][1];
                            }
/*
 * at this point, the first two equations read
 *              diag[0]*x0 +    u[0]*x1 + entry3*x2 = r[1]
 *       -ratio*diag[0]*x0 + diag[1]*x1 +   u[1]*x2 = 0
 * set r[0] = r[1] and eliminate x1 from the second equation
 */
                            s[3][0] = s[3][1];
           
                            s[1][1] += ratio*s[2][0];
                            s[2][1] += ratio*entry3;
                            s[3][1] = ratio*s[3][1];
                        } else {
/*
 * the i-th equation enforces continuity of the first derivative
 * across tau[i]; it reads
 *         -ratio*diag[i-1]*x_{i-1} + diag[i]*x_i + u[i]*x_{i+1} = r[i]
 * eliminate x_{i-1} from this equation
 */
                            s[1][i] += ratio*s[2][i - 1];
                            s[3][i] += ratio*s[3][i - 1];
                        }
/*
 * Set up the part of the next equation which depends on the i-th interval.
 */
                        if(z_half < 0) {
                            ratio = -s[5][i]*s[0][i]/s[1][i];
                            s[1][i + 1] = s[0][i]/3;
                        } else if(z_half == 0) {
                            ratio = -(s[0][i]/6)/s[1][i];
                            s[1][i + 1] = s[0][i]/3;
                        } else {
                            ratio = -(s[0][i]/6)/s[1][i];
                            s[1][i + 1] = s[0][i]*((1 - zeta*alpha)*factor/2 - s[5][i]);
                        }
                    }
    
                    s[4][ntau - 2] = 0.5;
/*
 * last two equations, which enforce continuity of third derivative and
 * of first derivative across tau[ntau - 2]
 */
                    double const entry_ = ratio*s[2][ntau - 3] + s[1][ntau - 2] + s[0][ntau - 2]/3;
                    s[1][ntau - 1] = s[0][ntau - 2]/6;
    s[3][ntau - 1] = ratio*s[3][ntau - 3] + s[3][ntau - 2];
    if (z_half < 0) {
       ratio = s[0][ntau - 2]*6*s[5][ntau - 3]*alpha/s[1][ntau - 3];
       s[1][ntau - 2] = ratio*s[2][ntau - 3] + s[0][ntau - 2] + s[0][ntau - 3];
       s[2][ntau - 2] = -s[0][ntau - 3];
    } else if (z_half == 0) {
       ratio = s[0][ntau - 2]/s[1][ntau - 3];
       s[1][ntau - 2] = ratio*s[2][ntau - 3] + s[0][ntau - 2] + s[0][ntau - 3];
       s[2][ntau - 2] = -s[0][ntau - 3];
    } else {
        const double factr2 = onemzt*(alpha*(onemzt*onemzt - 1) + 1)/(alpha*(onemzt*onemzt*onemzt - 1) + 1);
        ratio = factr2*s[0][ntau - 2]/s[1][ntau - 3];
        s[1][ntau - 2] = ratio*s[2][ntau - 3] + factr2*s[0][ntau - 3] + s[0][ntau - 2];
        s[2][ntau - 2] = -factr2*s[0][ntau - 3];
    }
/*
 * at this point, the last two equations read
 *             diag[i]*x_i +      u[i]*x_{i+1} = r[i]
 *      -ratio*diag[i]*x_i + diag[i+1]*x_{i+1} = r[i+1]
 *     eliminate x_i from last equation
 */
    s[3][ntau - 2] = ratio*s[3][ntau - 3];
    ratio = -entry_/s[1][ntau - 2];
    s[1][ntau - 1] += ratio*s[2][ntau - 2];
    s[3][ntau - 1] += ratio*s[3][ntau - 2];

/*
 * back substitution
 */
    s[3][ntau - 1] /= s[1][ntau - 1];
    for (int i = ntau - 2; i > 0; --i) {
        s[3][i] = (s[3][i] - s[2][i]*s[3][i + 1])/s[1][i];
    }

    s[3][0] = (s[3][0] - s[2][0]*s[3][1] - entry3*s[3][2])/s[1][0];
/*
 * construct polynomial pieces; first allocate space for the coefficients
 */
#if 1
/*
 * Start by counting the knots
 */
    {
        int const nknot0 = nknot;
        int nknot = ntau;
    
        for (int i = 0; i < ntau - 1; ++i) {
            double const z = s[4][i];
            if((z < 0.5 && z != 0.0) || (z > 0.5 && (1 - z) != 0.0)) {
                nknot++;
            }
        }
        assert (nknot == nknot0);
    }
#endif
    _allocateSpline(nknot);

    _knots[0] = tau[0];
    int j = 0;
    for (int i = 0; i < ntau - 1; ++i) {
        _coeffs[0][j] = gtau[i];
        _coeffs[2][j] = s[3][i];
        double const divdif = (gtau[i + 1] - gtau[i])/s[0][i];
        double z = s[4][i];
        double const z_half = z - 0.5;
        if (z_half < 0) {
           if (z == 0) {
              _coeffs[1][j] = divdif;
              _coeffs[2][j] = 0;
              _coeffs[3][j] = 0;
           } else {
              zeta = gamma*z;
              onemzt = 1 - zeta;
              double const c = s[3][i + 1]/6;
              double const d = s[3][i]*s[5][i];
              j++;
              
              double const del = zeta*s[0][i];
              _knots[j] = tau[i] + del;
              zt2 = zeta*zeta;
              double temp = onemg3/onemzt;
              alpha = (temp < 1 ? temp : 1);
              factor = onemzt*onemzt*alpha;
              temp = s[0][i];
              _coeffs[0][j] = gtau[i] + divdif*del + temp*temp*(d*onemzt*(factor - 1) + c*zeta*(zt2 - 1));
              _coeffs[1][j] = divdif + s[0][i]*(d*(1 - 3*factor) + c*(3*zt2 - 1));
              _coeffs[2][j] = (d*alpha*onemzt + c*zeta)*6;
              _coeffs[3][j] = (c - d*alpha)*6/s[0][i];
              if(del*zt2 == 0) {
                 _coeffs[1][j - 1] = 0; /* would be NaN in an */
                 _coeffs[3][j - 1] = 0; /*              0-length interval */
              } else {
                 _coeffs[3][j - 1] = _coeffs[3][j] - d*6*(1 - alpha)/(del*zt2);
                 _coeffs[1][j - 1] = _coeffs[1][j] -
                                   del*(_coeffs[2][j] - del/2*_coeffs[3][j - 1]);
              }
           }
        } else if (z_half == 0) {
           _coeffs[1][j] = divdif - s[0][i]*(s[3][i]*2 + s[3][i + 1])/6;
           _coeffs[3][j] = (s[3][i + 1] - s[3][i])/s[0][i];
        } else {
           onemzt = gamma*(1 - z);
           if (onemzt == 0) {
              _coeffs[1][j] = divdif;
              _coeffs[2][j] = 0;
              _coeffs[3][j] = 0;
           } else {
              zeta = 1 - onemzt;
              double const temp = onemg3/zeta;
              alpha = (temp < 1 ? temp : 1);
              double const c = s[3][i + 1]*s[5][i];
              double const d = s[3][i]/6;
              double const del = zeta*s[0][i];
              _knots[j + 1] = tau[i] + del;
              _coeffs[1][j] = divdif - s[0][i]*(2*d + c);
              _coeffs[3][j] = (c*alpha - d)*6/s[0][i];
              j++;

              _coeffs[3][j] = _coeffs[3][j - 1] +
                              (1 - alpha)*6*c/(s[0][i]*(onemzt*onemzt*onemzt));
              _coeffs[2][j] = _coeffs[2][j - 1] + del*_coeffs[3][j - 1];
              _coeffs[1][j] = _coeffs[1][j - 1] +
                               del*(_coeffs[2][j - 1] + del/2*_coeffs[3][j - 1]);
              _coeffs[0][j] = _coeffs[0][j - 1] +
                del*(_coeffs[1][j - 1] +
                            del/2*(_coeffs[2][j - 1] + del/3*_coeffs[3][j - 1]));
           }
        }

        j++;
        _knots[j] = tau[i + 1];
    }
/*
 * If there are discontinuities some of the knots may be at the same
 * position; in this case we generated some NaNs above. As they only
 * occur for 0-length segments, it's safe to replace them by 0s
 *
 * Due to the not-a-knot condition, the last set of coefficients isn't
 * needed (the last-but-one is equivalent), but it makes the book-keeping
 * easier if we _do_ generate them
 */
    double const del = tau[ntau - 1] - _knots[nknot - 2];
    
    _coeffs[0][nknot - 1] = _coeffs[0][nknot - 2] +
      del*(_coeffs[1][nknot - 2] + del*(_coeffs[2][nknot - 2]/2 +
                                       del*_coeffs[3][nknot - 2]/6));
    _coeffs[1][nknot - 1] = _coeffs[1][nknot - 2] +
      del*(_coeffs[2][nknot - 2] + del*_coeffs[3][nknot - 2]/2);
    _coeffs[2][nknot - 1] = _coeffs[2][nknot - 2] + del*_coeffs[3][nknot - 2];
    _coeffs[3][nknot - 1] = _coeffs[3][nknot - 2];

    assert (j + 1 == nknot);
}

void lsst::afw::math::detail::TautSpline::calculateTautSplineEvenOdd ( std::vector< double > const &  x, std::vector< double > const &  y, double const  gamma0, bool  even  )

private

Definition at line 1199 of file Spline.cc.

{
    const double *tau = &_tau[0];
    const double *gtau = &_gtau[0];
    int const ntau = _tau.size();       // size of tau and gtau, must be >= 2
    std::vector<double> x, y;           // tau and gtau, extended to -ve tau

    if(tau[0] == 0.0f) {
        int const np = 2*ntau - 1;
        x.resize(np);
        y.resize(np);

        x[ntau - 1] =  tau[0]; y[ntau - 1] = gtau[0];
        for (int i = 1; i != ntau; ++i) {
            if (even) {
                x[ntau - 1 + i] =  tau[i]; y[ntau - 1 + i] =  gtau[i];
                x[ntau - 1 - i] = -tau[i]; y[ntau - 1 - i] =  gtau[i];
            } else {
                x[ntau - 1 + i] =  tau[i]; y[ntau - 1 + i] =  gtau[i];
                x[ntau - 1 - i] = -tau[i]; y[ntau - 1 - i] = -gtau[i];
            }
        }
    } else {
        int const np = 2*ntau;
        x.resize(np);
        y.resize(np);

        for (int i = 0; i != ntau; ++i) {
            if (even) {
                x[ntau + i] =      tau[i]; y[ntau + i] =      gtau[i];
                x[ntau - 1 - i] = -tau[i]; y[ntau - 1 - i] =  gtau[i];
            } else {
                x[ntau + i] =      tau[i]; y[ntau + i] =      gtau[i];
                x[ntau - 1 - i] = -tau[i]; y[ntau - 1 - i] = -gtau[i];
            }
        }
    }

    TautSpline sp(x, y, gamma);         // fit a taut spline to x, y
/*
 * Now repackage that spline to reflect the original points
 */
    int ii;
    for (ii = sp._knots.size() - 1; ii >= 0; --ii) {
        if(sp._knots[ii] < 0.0f) {
            break;
        }
    }
    int const i0 = ii + 1;
    int const nknot = sp._knots.size() - i0;
   
    _allocateSpline(nknot);

    for (int i = i0; i != static_cast<int>(sp._knots.size()); ++i) {
        _knots[i - i0] = sp._knots[i];
        for (int j = 0; j != 4; ++j) {
            _coeffs[j][i - i0] = sp._coeffs[j][i];
        }
    }
}

---

The documentation for this class was generated from the following files:

• /home/lsstsw/stack/Linux64/afw/11.0-2-g04d2804/include/lsst/afw/math/detail/Spline.h
• /home/lsstsw/stack/Linux64/afw/11.0-2-g04d2804/src/math/Spline.cc