How I Stopped Worrying And Learned To Love Mathematica

A common tool in computer graphics are Catmull-Rom splines¹², which allow interpolating between a sequence of points. I used them in my dissertation, in which I created a handwriting recognition algorithm: For lines/strokes drawn on an iPad, I sampled the pen position and connected the points via these splines.

Since I was only interested in a qualitative description, I didn’t spend too much thought on their properties. When Freya Holmér asked on Twitter whether they are continuously differentiable, I was a bit disappointed in myself for not knowing. So, I sat down to find out.

As we’ll see below, working with general Catmull-Rom splines is a bit tedious. I started to compute their derivatives in my little notebook and soon got annoyed by it. I tried to bully Wolfram Alpha into giving me the answer, but the request was a bit too long and complicated for it. I needed the full Mathematica power.

Luckily, my university has a Mathematica licence for all their students and employees. Since I didn’t want to drive to campus just for that, I tried the remote desktop access. It was quite the hassle, but eventually, I made it work. Here’s what I found out… But first, the fundamentals on Catmull-Rom splines.³

The bare minimum on Catmull-Rom splines

Given a sequence of points \(P_0,P_1,...,P_n\in\mathbb{R}^2\) we might want find a continuous function \(\gamma:[0;1]\to\mathbb{R^2}\) that interpolates between the points. I.e., for some values

we want

to be exactly the points we are looking at. There are infinitely many ways to do so. For various reasons,⁴ a good way to do this is to find simple functions that connect two consecutive points and then stitch them together. Cubic polynomials of the form

are a good candidate for that.⁵ The whole function \(\gamma\) built from such cubic polynomials is then called a cubic spline.⁶

Let’s look at the four consecutive \(P_0, P_1, P_2, P_3\) somewhere within our actual (probably longer) sequence. First, we want that \(Z(b)=P_1\) and \(Z(c) = P_2\) for it to be an interpolation of the middle bit. Then, we want \(P_0, P_3\) to somehow define the exact shape of this polynomial. There are still many ways to do this,⁷⁸ and here is one created by Edwin Catmull and Raphael Rom⁹ in a form presented by Phillip Barry and Ronald Goldman¹⁰:

For an \(\alpha\in[0;1]\) we recursively define four knot weights \(t_0, t_1, t_2, t_3\) via

I.e., we take the Euclidean distances of consecutive points \(\lVert P_{i+1}-P_i\rVert_2\) and add them up. But we take them to the power of \(\alpha\) to be able to control the shape of the resulting curve.¹¹ We use these knot weights to interpolate between various auxiliary points. To do so, let’s abbreviate¹²

for \(s\in\mathbb{R}\) and \(A,B\in\mathbb{R^2}\). Then, for \(t\in[t_1, t_2]\) we define…

There’s much more one could say about this and why this leads to nice and useful splines. But for the rest of this post, we’re only concerned with their derivatives with respect to \(t\): When the whole spline is used for some sort of road, rollercoaster track, etc. the transitions from one part to the next should be smooth. I.e., the first derivative of the cubic polynomials used should be the same at the points \(P_i\) in which they meet. This way, there will be no sudden velocity changes when you follow the curve in a vehicle. Moreover, equal higher derivatives are also often desirable because of similar physical considerations. So, let’s see how they look like.

Derivatives of Catmull-Rom splines

First, we rename the points \(P_0,P_1,P_2,P_3\) as \(P,Q,R,S\) and the weights \(t_0,t_1,t_2,t_3\) as \(a,b,c,d\) to make things slightly easier to read. In particular, I used these letters in my Mathematica code:

    In[1]: u[p_,q_,r_,s_,a_,b_,c_,d_,t_] := (b-t)/(b-a) * p + (t-a)/(b-a) * q

    In[2]: v[p_,q_,r_,s_,a_,b_,c_,d_,t_] := (c-t)/(c-b) * q + (t-b)/(c-b) * r

    In[3]: w[p_,q_,r_,s_,a_,b_,c_,d_,t_] := (d-t)/(d-c) * r + (t-c)/(d-c) * s

    In[4]: x[p_,q_,r_,s_,a_,b_,c_,d_,t_] := (c-t)/(c-a) * u[p,q,r,s,a,b,c,d,t] + (t-a)/(c-a) * v[p,q,r,s,a,b,c,d,t]

    In[5]: y[p_,q_,r_,s_,a_,b_,c_,d_,t_] := (d-t)/(d-b) * v[p,q,r,s,a,b,c,d,t] + (t-b)/(d-b) * w[p,q,r,s,a,b,c,d,t]

    In[6]: z[p_,q_,r_,s_,a_,b_,c_,d_,t_] := (c-t)/(c-b) * x[p,q,r,s,a,b,c,d,t] + (t-b)/(c-b) * y[p,q,r,s,a,b,c,d,t]

This and Mathematica’s TexForm function, which spits out LaTeX code, allow us to find a singular term for the cubic polynomial \(Z\) in question:

Feel free to open it in a new tab… Nice, isn’t it?

Since everything is sorted with respect to powers of \(t\), it’s straightforward to compute the derivative. In general, the result won’t look much different, though. To get closer to our goal of understanding how two consecutive cubic polynomials fit together, we have to analyse them at their endpoints.

Let’s assume this \(Z\) introduced above is the first of two such parts. Since it interpolates between \(Q\) and \(R\), we want to look at its derivative at \(t=c\). (Because \(Z(c) = R\), which Mathematica verifies for us, just to make sure.) Moreover, we can look at the knot weights a bit closer again: The exact form obtained via distances of points isn’t that important here. But it’s crucial that they are cumulative sums. I.e., we can find \(e,f,g>0\) such that

With this, we can feed the line

    In[7]: Simplify[Expand[D[z[P, Q, R, S, 0, e, e+f, e+f+g, t], t] /. t -> e + f]]

into Mathematica to obtain

The next cubic polynomial in the spline uses the points \(Q,R,S\) and a new point \(T\). The knot weights are also shifted to include a new value \(h>0\). Explicitly, they now are

Moreover, we now want to analyse the polynomial at the start, i.e., at \(t=f\). So, with

    In[8]: Simplify[Expand[D[z[Q, R, S, T, 0, f, f+g, f+g+h, t], t] /. t -> f]]

we get

And it’s the same as the first one!

This means that two cubic polynomials built this way will fit together smoothly. In other words, Catmull-Rom splines are always once continuously differentiable. This was to be expected just from looking at them (e.g. the gif at the top of this post). But it’s nice to see this confirmed. In particular, we only used the fact, that the knot weights are cumulative sums of non-negative numbers. Their exact form is irrelevant.

The next step would be to look at the second derivatives, and without further ado we get:

and

If you can’t see it directly: In the first one we have \(P\) and in the second one we don’t. A bit crude, but basically enough to say they aren’t the same in general. So, Catmull-Rom splines aren’t two times continuously differentiable.

Final note

There’s much more to explore here.

First, the form of the first derivate

is quite interesting. There’s a symmetry in the summands which can be interpreted geometrically. Basically, an even mix between where the curve is comming from and where it is going.

Second, one could surely endeavour to find all knot weights and \(\alpha\)-values for which the second derivates do coincide. From a quick look it seems that \(f\) or \(g\) being \(0\) might be a necessary condition; which would mean two consecutive points to coincide. But needs more time…

This post is long enough, though, so let’s end it here.

Final final note

Was this a weird first blog post? Should I have written a more detailed primer on Catmull-Rom splines first? Should I have chosen a different title when the post isn’t really about Mathematica?

Will I write more coherent stuff in the future?

I mean… I’ll try!