Smooth Paths Using Catmull-Rom Splines

Let’s assume that we have a set of points and we want to draw a path that goes through them. The simplest and easiest way to do this is to connect the points with straight lines as shown in Figure 1.

A path using straight lines.

However, straight lines would cause sharp corners to the path and in many cases a smoother path would be preferred. Figure 2 shows a smoother path that goes through the same points that are in Figure 1.

A smooth path.

In this post I’m going to focus on Catmull-Rom splines which are commonly used in computer graphics to create smooth curves. For example, the path in Figure 2 uses them.

Catmull-Rom splines

Catmull-Rom splines are piecewise-defined polynomial functions. A single spline segment is defined by four control points $\boldsymbol{p}_0, \dots, \boldsymbol{p}_3$ but the actual curve is drawn only between points $\boldsymbol{p}_1$ and $\boldsymbol{p}_2$ as is illustrated in Figure 3. However, it is easy to chain these segments together.

One segment of Catmull-Rom spline.

If we want to draw a curve that goes through $k$ points, we need $k+2$ control points because the curve is not drawn through the first and the last ones. These two additional points can be selected arbitrarily, but they affect the shape of the curve. Now a segment between points $\boldsymbol{p}_n$ and $\boldsymbol{p}_{n+1}$ is calculated using points $\boldsymbol{p}_{n-1}$, $\boldsymbol{p}_n$, $\boldsymbol{p}_{n+1}$, and $\boldsymbol{p}_{n+2}$, where $1 \leq n \leq k - 1$. When these segments are combined together, they form a continuous curve, which passes through all points between $\boldsymbol{p}_1$ and $\boldsymbol{p}_{k}$. Figure 2 shows an example of this kind of curve.

There are some parameters that can be used to control the shape of the spline. Catmull-Rom splines have three common variants: uniform, centripetal, and chordal. Differences have been studied for example here. Additionally, it is possible to use a tension parameter that defines how “tight” the spline is. When tension is 0, the result looks like the curve in Figure 2, and when tension is 1, the result is straight lines as in Figure 1. The following interactive example demonstrates Catmull-Rom splines and these parameters.

You can add points by clicking the canvas and move existing points by dragging them. The control panel allows you to change the tension of the curve and select the spline variant that is drawn.

Calculating spline segments

Let’s start with following equations that are described in a Wikipedia article about centripetal Catmull-Rom splines:

$$\boldsymbol{q}(t) = \frac{t_2 - t}{t_2 - t_1}\boldsymbol{B}_1 + \frac{t - t_1}{t_2 - t_1}\boldsymbol{B}_2,$$

where

$$\begin{align} \boldsymbol{B}_1 &= \frac{t_2 - t}{t_2 - t_0}\boldsymbol{A}_1 + \frac{t - t_0}{t_2 - t_0}\boldsymbol{A}_2, \\ \boldsymbol{B}_2 &= \frac{t_3 - t}{t_3 - t_1}\boldsymbol{A}_2 + \frac{t - t_1}{t_3 - t_1}\boldsymbol{A}_3, \\ \boldsymbol{A}_1 &= \frac{t_1 - t}{t_1 - t_0}\boldsymbol{p}_0 + \frac{t - t_0}{t_1 - t_0}\boldsymbol{p}_1, \\ \boldsymbol{A}_2 &= \frac{t_2 - t}{t_2 - t_1}\boldsymbol{p}_1 + \frac{t - t_1}{t_2 - t_1}\boldsymbol{p}_2, \\ \boldsymbol{A}_3 &= \frac{t_3 - t}{t_3 - t_2}\boldsymbol{p}_2 + \frac{t - t_2}{t_3 - t_2}\boldsymbol{p}_3, \end{align}$$

and where

$$t_i = t_{i-1} + \|\boldsymbol{p}_i - \boldsymbol{p}_{i-1}\|^\alpha,$$

and $t_0 = 0$, $i = 1, 2, 3$, and $\alpha \in [0, 1]$.

The actual segment we are interested in is between $t_1$ and $t_2$ i.e. $\boldsymbol{q}(t_1) = \boldsymbol{p}_1$ and $\boldsymbol{q}(t_2) = \boldsymbol{p}_2$. If $\alpha = 0$, the resulting curve $\boldsymbol{q}$ is a uniform Catmull-Rom spline. When $\alpha = 0.5$, the curve is a centripetal variant and when $\alpha = 1$, the result is a chordal variant.

Calculating the curve with previous equation can be quite inconvenient. Often it would be better to use precalculated constants $\boldsymbol{a}$, $\boldsymbol{b}$, $\boldsymbol{c}$, and $\boldsymbol{d}$, and represent the curve segment between $\boldsymbol{p}_1$ and $\boldsymbol{p}_2$ as

$$\boldsymbol{p}(t) = \boldsymbol{a}t^3 + \boldsymbol{b}t^2 + \boldsymbol{c}t + \boldsymbol{d},$$

where $t \in [0, 1]$, $\boldsymbol{p}(0) = \boldsymbol{p}_1$, and $\boldsymbol{p}(1) = \boldsymbol{p}_2$. Other way to define this is

$$\boldsymbol{p}(t) = \boldsymbol{q}(t_1 + t(t_2 - t_1)).$$

We can get a relationship between the tangents of $\boldsymbol{p}$ and $\boldsymbol{q}$ by differentiating these with respect to $t$ as follows:

$$\boldsymbol{p}'(t) = 3\boldsymbol{a}t^2 + 2\boldsymbol{b}t + \boldsymbol{c} = (t_2 - t_1)\boldsymbol{q}'(t_1 + t(t_2 - t_1)).$$

Now we have

$$\begin{align} \boldsymbol{p}(0) &= \boldsymbol{q}(t_1) &&= \boldsymbol{p}_1 = \boldsymbol{d}, \\ \boldsymbol{p}(1) &= \boldsymbol{q}(t_2) &&= \boldsymbol{p}_2 = \boldsymbol{a} + \boldsymbol{b} + \boldsymbol{c} + \boldsymbol{d}, \\ \boldsymbol{p}'(0) &= (t_2 - t_1)\boldsymbol{q}'(t_1) &&= \boldsymbol{m}_1 = \boldsymbol{c}, \\ \boldsymbol{p}'(1) &= (t_2 - t_1)\boldsymbol{q}'(t_2) &&= \boldsymbol{m}_2 = 3\boldsymbol{a} + 2\boldsymbol{b} + c, \end{align}$$

where symbols $\boldsymbol{m}_1$, and $\boldsymbol{m}_2$ are used to represent tangents at the starting point $\boldsymbol{p}_1$ and at the ending point $\boldsymbol{p}_2$ respectively. By solving $\boldsymbol{a}$, $\boldsymbol{b}$, $\boldsymbol{c}$, and $\boldsymbol{d}$ from these four equations we get

$$\begin{align} \boldsymbol{a} &= 2\boldsymbol{p}_1 - 2\boldsymbol{p}_2 + \boldsymbol{m}_1 + \boldsymbol{m}_2, \\ \boldsymbol{b} &= -3\boldsymbol{p}_1 + 3\boldsymbol{p}_2 - 2\boldsymbol{m}_1 - \boldsymbol{m}_2, \\ \boldsymbol{c} &= \boldsymbol{m}_1, \\ \boldsymbol{d} &= \boldsymbol{p}_1. \end{align}$$

Now we have to determine $\boldsymbol{m}_1$ and $\boldsymbol{m}_2$ and for that we need tangents $\boldsymbol{q}'(t_1)$ and $\boldsymbol{q}'(t_2)$. Calculating the derivative of $\boldsymbol{q}$ is quite cumbersome, but luckily we can use Mathematica to do the work for us. We get

$$\begin{align} \boldsymbol{q}'(t_1) &= \frac{\boldsymbol{p}_1 - \boldsymbol{p}_0}{t_1 - t_0} - \frac{\boldsymbol{p}_2 - \boldsymbol{p}_0}{t_2 - t_0} + \frac{\boldsymbol{p}_2 - \boldsymbol{p}_1}{t_2 - t_1}, \\ \boldsymbol{q}'(t_2) &= \frac{\boldsymbol{p}_2 - \boldsymbol{p}_1}{t_2 - t_1} - \frac{\boldsymbol{p}_3 - \boldsymbol{p}_1}{t_3 - t_1} + \frac{\boldsymbol{p}_3 - \boldsymbol{p}_2}{t_3 - t_2}, \end{align}$$

and thus

$$\begin{align} \boldsymbol{m}_1 = (t_2 - t_1)\left(\frac{\boldsymbol{p}_1 - \boldsymbol{p}_0}{t_1 - t_0} - \frac{\boldsymbol{p}_2 - \boldsymbol{p}_0}{t_2 - t_0} + \frac{\boldsymbol{p}_2 - \boldsymbol{p}_1}{t_2 - t_1}\right), \\ \boldsymbol{m}_2 = (t_2 - t_1)\left(\frac{\boldsymbol{p}_2 - \boldsymbol{p}_1}{t_2 - t_1} - \frac{\boldsymbol{p}_3 - \boldsymbol{p}_1}{t_3 - t_1} + \frac{\boldsymbol{p}_3 - \boldsymbol{p}_2}{t_3 - t_2}\right). \end{align}$$

If we want to take tension $\tau$ into account, we must also multiply both $\boldsymbol{m}_1$ and $\boldsymbol{m}_2$ with $1 - \tau$.

It is now easy to precalculate coefficients $\boldsymbol{a}$, $\boldsymbol{b}$, $\boldsymbol{c}$, and $\boldsymbol{d}$ for each segment and use equation $\boldsymbol{p}(t) = \boldsymbol{a}t^3 + \boldsymbol{b}t^2 + \boldsymbol{c}t + \boldsymbol{d}$ to quickly interpolate that segment.

C++ implementation

The code is written in C++. It uses GLM library that provides the basic linear algebra functionality. Overall, I’ll try to keep the code syntax in this blog quite simple so that it is easy to understand and to convert to other languages.

In its simplest form, we can define a struct of the spline segment as follows:

struct Segment
{
    vec2 a;
    vec2 b;
    vec2 c;
    vec2 d;
};

Note that we are using two-dimensional splines here. If you want three-dimensional splines, just replace all occurances of vec2 with vec3.

As said in the previous chapter, we need parameters $\alpha$ and $\tau$ to calculate the spline. In following code we use variable alpha for $\alpha$ and tension for $\tau$. A good value for alpha is 0.5 which gives us a centripetal Catmull-Rom spline, and for tension a value 0 is a good choice. These values can range from 0 to 1.

If we now have points p0, p1, p2, and p3 for each segment, we can calculate coefficients for segments with equations from previous chapter:

float t0 = 0.0f;
float t1 = t0 + pow(distance(p0, p1), alpha);
float t2 = t1 + pow(distance(p1, p2), alpha);
float t3 = t2 + pow(distance(p2, p3), alpha);

vec2 m1 = (1.0f - tension) * (t2 - t1) *
    ((p1 - p0) / (t1 - t0) - (p2 - p0) / (t2 - t0) + (p2 - p1) / (t2 - t1));
vec2 m2 = (1.0f - tension) * (t2 - t1) *
    ((p2 - p1) / (t2 - t1) - (p3 - p1) / (t3 - t1) + (p3 - p2) / (t3 - t2));

Segment segment;
segment.a = 2.0f * (p1 - p2) + m1 + m2;
segment.b = -3.0f * (p1 - p2) - m1 - m1 - m2;
segment.c = m1;
segment.d = p1;

We can get the same result slightly more efficiently by simplifying the equations and using the following code to calculate m1 and m2:

float t01 = pow(distance(p0, p1), alpha);
float t12 = pow(distance(p1, p2), alpha);
float t23 = pow(distance(p2, p3), alpha);

vec2 m1 = (1.0f - tension) *
    (p2 - p1 + t12 * ((p1 - p0) / t01 - (p2 - p0) / (t01 + t12)));
vec2 m2 = (1.0f - tension) *
    (p2 - p1 + t12 * ((p3 - p2) / t23 - (p3 - p1) / (t12 + t23)));

These coefficients can be precaluclated for all spline segments assuming that the parameters do not change afterwards. It is now simple to retrieve a point in a segment segment by using the precalculated values:

vec2 point = segment.a * t * t * t +
             segment.b * t * t +
             segment.c * t +
             segment.d;

And t is of course a value between 0 and 1. With a value 0, the result is the starting point of the spline segment, and with a value 1, the result is the ending point.