Natural cubic spline

Piecewise spline interpolation fits cubic polynomials through a set of points. In contrast to utilizing a polynomial of a higher degree, this results in a smooth interpolation that stays much closer to the target points. The interpolation is based on a t value. The t value has to be monotonically increasing. This in turn means that a regular spline can only ever go into one direction. However, parametric splines overcome this limitation.

The math

The construction of a natural cubic spline is described in “Numerical Analysis” by Richard L. Burden and J. Douglas Faires. In my implementation the t value increases by one for every segment.

void ComputeCubicSpline(double* values, int valueCount)
{
    double* x = new double[valueCount];
    double* a = new double[valueCount];
    for (int i = 0; i < valueCount; i++)
    {
		// segment t values: [0,1], [1,2] This is asumed when finding the segment to a t value
        x[i] = i; 
        a[i] = values[i];
    }

    if (valueCount == 0)
    {
        printf("Error: Trying to compute spline without points\n");
        return;
    }

    // number of segments is equal to number of points - 1
    int segmentCount = valueCount - 1;

    // Implementation of the books algorithm
    // allocate memory
    double* b = new double[segmentCount + 1];
    double* c = new double[segmentCount + 1];
    double* d = new double[segmentCount + 1];

    // step 1
    double* h = new double[segmentCount];
    for (int i = 0; i <= segmentCount - 1; i++)
    {
        h[i] = x[i + 1] - x[i];
    }

    // step 2
    double* alpha = new double[segmentCount];
    for (int i = 1; i <= segmentCount - 1; i++)
    {
        alpha[i] = (3 / h[i]) * (a[i + 1] - a[i]) - (3 / h[i - 1]) * (a[i] - a[i - 1]);
    }

    // step 3
    // step 3,4,5 and part of 6 solve tridiagonal system
    double* l = new double[segmentCount + 1];
    double* u = new double[segmentCount + 1];
    double* z = new double[segmentCount + 1];
    l[0] = 1;
    u[0] = 0;
    z[0] = 0;

    // step 4
    for (int i = 1; i <= segmentCount - 1; i++)
    {
        l[i] = 2 * (x[i + 1] - x[i - 1]) - h[i - 1] * u[i - 1];
        u[i] = h[i] / l[i];
        z[i] = (alpha[i] - h[i - 1] * z[i - 1]) / l[i];
    }

    delete[] x;
    delete[] alpha;

    // step 5
    l[segmentCount] = 1;
    z[segmentCount] = 0;
    c[segmentCount] = 0;

    delete[] l;

    // step 6
    for (int j = segmentCount - 1; j >= 0; j--)
    {
        c[j] = z[j] - u[j] * c[j + 1];
        b[j] = (a[j + 1] - a[j]) / h[j] - h[j] * (c[j + 1] + 2 * c[j]) / 3;
        d[j] = (c[j + 1] - c[j]) / (3 * h[j]);
    }

    delete[] h;
    delete[] u;
    delete[] z;

    for (int i = 0; i < segmentCount; i++)
    {
		// store polynomials
        polynomial_.push_back(Polynomial(d[i], c[i], b[i], a[i]));
    }

    delete[] a;
    delete[] b;
    delete[] c;
    delete[] d;
}

To get a position on the spline we just enter a t value into the corresponding polynomial. Since we know that t is increased by one for every segment, we can directly use t to lookup the corresponding segment/polynomial.

Parametric spline

Parametric splines are based on a cubic spline for each dimension. Therefore, a 2D spline consists of two cubic splines and a 3D spline of three.

void ComputeParametricSpline(double* x, double* y, double* z, int valueCount)
{
    splineX_ = CubicSpline();
    splineY_ = CubicSpline();
    splineZ_ = CubicSpline();

    splineX_.Compute(x, valueCount);
    splineY_.Compute(y, valueCount);
    splineZ_.Compute(z, valueCount);
}

To get the position on a three dimensional parametric spline, we just pass our t value to each of our splines:

Vector3 eval(float t)
{
    return Vector3(splineX_.eval(t), splineY_.eval(t), splineZ_.eval(t));
}

Since we now have a cubic spline for each dimension, we can move freely into any direction and also intersect ourselves.

Demo inside my 2D engine

This demo constructs a 2D parametric spline from points added by clicking. Furthermore, it animates a point on the spline by passing the current time as t value.

Prototype implementation

My parametric spline implementation is available on my GitHub.