# Cubic Spline Interpolation

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:

The code for the spline class is below:

# 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.

- Create new arrays $a, b, d$ of size $n$
- Create new arrays $ c, h, l, u, z$ of size $ n+1$
- Set $l_0 = 1, u_0 = z_0 = 0$
- Set $h_0 = x_1 - x_0$
- For $i = 1, \dots, n-1$
- Set $h_i = x_{i+1} - x_i$
- Set $l_i = 2 \left ( x_{i+1} - x_{i-1} \right ) - h_{i-1}u_{i-1}$
- Set $u_i = \frac{h_i}{l_i}$
- 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 )$
- Set $z_i = \frac{a_i - h_{i-1}z_{i-1}}{l_i}$

- Set $l_n = 1, z_n = c_n = 0$
- For $j = n-1, n-2, \dots, 0$
- Set $c_j = z_j - u_j c_{j+1}$
- Set $b_j = \frac{y_{j+1} - y_j}{h_j} - \frac{h_j \left( c_{j+1} + 2c_j \right)}{3}$
- Set $d_j = \frac{c_{j+1} - c_j}{3h_j}$

- For $i = 0, \dots, n-1$
- Set $S_{i,\left(x_i,a,b,c,d\right)} = \left( x_i, y_i, b_i, c_i, d_i \right)$

12 comments:

Linzy Cumbia at 2012-02-23 23:07:18 -0500

How would I get results out of this to use them to draw a graph on a UIView in Objective-C?

Il Jae Lee at 2013-01-28 08:58:00 -0500

Thank you so much for posting such code. This will be so much helpful.

Giovanni at 2013-06-11 04:35:25 -0400

Thanks for posting this! Finally after more than one month I've finally been able to implement the CubicSpline interpolation!!!

Guy at 2013-07-03 11:47:51 -0400

Thank you for the code ! It works great

khan at 2013-10-29 05:11:41 -0400

How one can use that class. Could any body help to provide the main code. regards

Ali Alsaqqa at 2014-07-11 00:57:52 -0400

Thanks. Is there any usage example?

Victor at 2014-08-26 20:24:26 -0400

Thanks Man! This is the best spline library I saw today and I saw quite a few. I'll definitely buy you a beer!

Mark Mitchell at 2015-03-14 23:56:07 -0400

Thanks. Good stuff. I'm trying out a modified version of your library that add support for calculating bezier control points. You can find it at github (Spline directory of Engauge6). Some parts may be useful to you, although I'm not recommending you remove the templates like I did for my very specific use case. BTW, any suggestions for delivering on the beer fee?

mviljamaa at 2015-12-08 02:09:42 -0500

Could you clarify how is this 2D? You have times (t-axis) and two points (x,y) in a vec2?

mviljamaa at 2015-12-08 02:19:46 -0500

Also, What does the interpolate function do?

mviljamaa at 2015-12-08 02:41:29 -0500

Also I have the following problem: Error: no suitable conversion function from Spline to "float" exists When trying to pass sp[x] to a function taking float.

mviljamaa at 2015-12-08 03:29:17 -0500

Is there some alternative data type that one could use instead of glm::vec2?