Cubic spline interpolation is a simple way to obtain a smooth curve from a set of discrete points (knots). It has both C1 (first derivative) and C2 (second derivative) continuity, enabling it to produce a continuous piecewise function given a set of data points.

From the algorithm detailed below I have implemented a clamped cubic spline class in C++. It is templated on the type of X and Y, allowing for use of scalar or vector types. It requires only that X and Y define the operators +, -, *, and /, and that Y have a constructor that takes a single scalar argument, initializing all elements to the supplied value.

Usage is very simple. For example, to interpolate a 2D location over time, try this:

#import "glm.hpp"
#import <vector>

std::vector<float> times;
std::vector<glm::vec2> points;
times.push_back(0);
points.push_back(glm::vec2(0));
times.push_back(.5f);
points.push_back(glm::vec2(-303, -572));
times.push_back(1);
points.push_back(glm::vec2(-250, -980));

/* Create the spline interpolating the position over time */
Spline<float, glm::vec2> sp(times, points);

/* Compute the interpolated position at time 0.75f */
glm::vec2 value(sp[.75f]);

The code for the spline class is below:

/* "THE BEER-WARE LICENSE" (Revision 42): Devin Lane wrote this file. As long as you retain 
 * this notice you can do whatever you want with this stuff. If we meet some day, and you
 * think this stuff is worth it, you can buy me a beer in return. */

#include <vector>
#include <iostream>

/** Templated on type of X, Y. X and Y must have operator +, -, *, /. Y must have defined
 * a constructor that takes a scalar. */
template <typename X, typename Y>
class Spline {
public:
    /** An empty, invalid spline */
    Spline() {}
    
    /** A spline with x and y values */
    Spline(const std::vector<X>& x, const std::vector<Y>& y) {
        if (x.size() != y.size()) {
            std::cerr << "X and Y must be the same size " << std::endl;
            return;
        }
        
        if (x.size() < 3) {
            std::cerr << "Must have at least three points for interpolation" << std::endl;
            return;
        }
        
        typedef typename std::vector<X>::difference_type size_type;
        
        size_type n = y.size() - 1;
        
        std::vector<Y> b(n), d(n), a(n), c(n+1), l(n+1), u(n+1), z(n+1);
        std::vector<X> h(n+1);

        l[0] = Y(1);
        u[0] = Y(0);
        z[0] = Y(0);
        h[0] = x[1] - x[0];
            
        for (size_type i = 1; i < n; i++) {
            h[i] = x[i+1] - x[i];
            l[i] = Y(2 * (x[i+1] - x[i-1])) - Y(h[i-1]) * u[i-1];
            u[i] = Y(h[i]) / l[i];
            a[i] = (Y(3) / Y(h[i])) * (y[i+1] - y[i]) - (Y(3) / Y(h[i-1])) * (y[i] - y[i-1]);
            z[i] = (a[i] - Y(h[i-1]) * z[i-1]) / l[i];
        }
            
        l[n] = Y(1);
        z[n] = c[n] = Y(0);
        
        for (size_type j = n-1; j >= 0; j--) {
            c[j] = z[j] - u[j] * c[j+1];
            b[j] = (y[j+1] - y[j]) / Y(h[j]) - (Y(h[j]) * (c[j+1] + Y(2) * c[j])) / Y(3);
            d[j] = (c[j+1] - c[j]) / Y(3 * h[j]);
        }
        
        for (size_type i = 0; i < n; i++) {
            mElements.push_back(Element(x[i], y[i], b[i], c[i], d[i]));
        }        
    }
    virtual ~Spline() {}
    
    Y operator[](const X& x) const {
        return interpolate(x);
    }
    
    Y interpolate(const X&x) const {
        if (mElements.size() == 0) return Y();
        
        typename std::vector<element_type>::const_iterator it;
        it = std::lower_bound(mElements.begin(), mElements.end(), element_type(x));
        if (it != mElements.begin()) {
            it--;
        }   
            
        return it->eval(x);
    }
    
    std::vector<Y> operator[](const std::vector<X>& xx) const {
        return interpolate(xx);
    }
    
    /* Evaluate at multiple locations, assuming xx is sorted ascending */
    std::vector<Y> interpolate(const std::vector<X>& xx) const {
        if (mElements.size() == 0) return std::vector<Y>(xx.size());
        
        typename std::vector<X>::const_iterator it;
        typename std::vector<element_type>::const_iterator it2;
        it2 = mElements.begin();
        std::vector<Y> ys;
        for (it = xx.begin(); it != xx.end(); it++) {
            it2 = std::lower_bound(it2, mElements.end(), element_type(*it));
            if (it2 != mElements.begin()) {
                it2--;
            }
                
            ys.push_back(it2->eval(*it));
        }

        return ys;
    }

protected:
    
    class Element {
    public:
        Element(X _x) : x(_x) {}
        Element(X _x, Y _a, Y _b, Y _c, Y _d)
        : x(_x), a(_a), b(_b), c(_c), d(_d) {}
        
        Y eval(const X& xx) const {
            X xix(xx - x);
            return a + b * xix + c * (xix * xix) + d * (xix * xix * xix);
        }
        
        bool operator<(const Element& e) const {
            return x < e.x;
        }
        bool operator<(const X& xx) const {
            return x < xx;
        }
        
        X x;
        Y a, b, c, d;
    };
            
    typedef Element element_type;
    std::vector<element_type> mElements;
};

Math

Spline interpolation is useful for creating “continuous” curves out of discrete data points. These curves are can have C1 (first derivative) and C2 (second derivative) continuity, allowing for a smooth join between segments of the curve. An algorithm for computing natural cubic splines was listed on Wikipedia, and have reproduced it here so as to have the full math -> code translation in one place.

Definition of a Cubic Spline

Given a list of data values

where $y_n$ is an n-tuple defined as

with corresponding location values $X$ defined as

we wish to compute a value of Y for any value in the range $[x_0, x_n]$. We will do this for some $x_i \leq x < x_j$ by evaluating a cubic spline $S_i(x)$.

A cubic spline is of the form:

where $(x_i, a, b, c, d)$ is a 5-tuple describing the parameters of $S_i(x)$. Given our set of data $Y$ and locations $X$, we wish to find $n$ polynomials $S_i (x)$ for $i = 0, \dots , n-1$ such that

We can compute the parameters of these $n$ polynomials using the algorithm below.

Computation of Cubic Spline Parameters

Input: list of points $Y$, with $\left| Y \right| =n+1$, with locations $X$

Output: list spline parameters composed of $n$ 5-tuples.

  1. Create new arrays $a, b, d$ of size $n$
  2. Create new arrays $ c, h, l, u, z$ of size $ n+1$
  3. Set $l_0 = 1, u_0 = z_0 = 0$
  4. Set $h_0 = x_1 - x_0$
  5. For $i = 1, \dots, n-1$
    1. Set $h_i = x_{i+1} - x_i$
    2. Set $l_i = 2 \left ( x_{i+1} - x_{i-1} \right ) - h_{i-1}u_{i-1}$
    3. Set $u_i = \frac{h_i}{l_i}$
    4. Set $a_i = \frac{3}{h_i} \left ( y_{i+1} - y_i \right ) - \frac{3}{h_{i-1}} \left ( y_i - y_{i-1} \right )$
    5. Set $z_i = \frac{a_i - h_{i-1}z_{i-1}}{l_i}$
  6. Set $l_n = 1, z_n = c_n = 0$
  7. For $j = n-1, n-2, \dots, 0$
    1. Set $c_j = z_j - u_j c_{j+1}$
    2. Set $b_j = \frac{y_{j+1} - y_j}{h_j} - \frac{h_j \left( c_{j+1} + 2c_j \right)}{3}$
    3. Set $d_j = \frac{c_{j+1} - c_j}{3h_j}$
  8. For $i = 0, \dots, n-1$
    1. Set $S_{i,\left(x_i,a,b,c,d\right)} = \left( x_i, y_i, b_i, c_i, d_i \right)$