Class TautSpline

Inheritance Relationships

Base Type

Class Documentation

class TautSpline : public lsst::afw::math::detail::Spline

Public Types

enum Symmetry

Values:

Unknown
Odd
Even

Public Functions

TautSpline(std::vector<double> const &x, std::vector<double> const &y, double const gamma = 0, Symmetry type = Unknown)

Construct cubic spline interpolant to given data.

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)

If gamma > 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 notes below. The interpolant gets rounded with increasing gamma. A value of 2.5 for gamma is typical.

  • in (3, 6), same as above, 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.

Note

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.

Note

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.

Note

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)

Parameters

  • x: points where function’s specified

  • y: values of function at tau[]

  • gamma: control extra knots. See main description for details.

  • type: specify the desired symmetry (e.g. Even)

Exceptions

  • pex::exceptions::InvalidParameterError: Thrown if x and y do not have the same length or do not have at least two points