How I Stopped Worrying And Learned To Love Mathematica
10th March 2022
A common tool in computer graphics are CatmullRom splines^{1}^{2}, 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 CatmullRom 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 CatmullRom splines.^{3}
The bare minimum on CatmullRom 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
\[0=t_0 <\cdots < t_n =1\]we want
\[\gamma(t_0) = P_0,\quad \gamma(t_1)=P_1,\quad ...,\quad \gamma(t_n)=P_n\]to be exactly the points we are looking at. There are infinitely many ways to do so. For various reasons,^{4} 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
\[Z:[b,c]\to\mathbb{R^2}, \quad t\mapsto At^3 + Bt^2 + Ct + D\]are a good candidate for that.^{5} The whole function \(\gamma\) built from such cubic polynomials is then called a cubic spline.^{6}
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,^{7}^{8} and here is one created by Edwin Catmull and Raphael Rom^{9} in a form presented by Phillip Barry and Ronald Goldman^{10}:
For an \(\alpha\in[0;1]\) we recursively define four knot weights \(t_0, t_1, t_2, t_3\) via
\[t_0 = 0,\] \[t_{i+1} = t_i + \lVert P_{i+1}P_i\rVert_2^\alpha\qquad \text{for } i=1,2,3.\]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.^{11} We use these knot weights to interpolate between various auxiliary points. To do so, let’s abbreviate^{12}
\[\mathrm{lerp}(A,B,s)\enspace :=\enspace s\cdot B + (1s)\cdot A\]for \(s\in\mathbb{R}\) and \(A,B\in\mathbb{R^2}\). Then, for \(t\in[t_1, t_2]\) we define…
\[U := \mathrm{lerp}\left(P_0, P_1, \frac{tt_0}{t_1t_0}\right),\] \[V := \mathrm{lerp}\left(P_1, P_2, \frac{tt_1}{t_2t_1}\right),\] \[W := \mathrm{lerp}\left(P_2, P_3, \frac{tt_2}{t_3t_2}\right),\] \[X := \mathrm{lerp}\left(U, V, \frac{tt_1}{t_3t_1}\right),\] \[Y := \mathrm{lerp}\left(V, W, \frac{tt_2}{t_4t_2}\right),\] \[Z := \mathrm{lerp}\left(X, Y, \frac{tt_2}{t_3t_2}\right).\]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 CatmullRom 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_] := (bt)/(ba) * p + (ta)/(ba) * q
In[2]: v[p_,q_,r_,s_,a_,b_,c_,d_,t_] := (ct)/(cb) * q + (tb)/(cb) * r
In[3]: w[p_,q_,r_,s_,a_,b_,c_,d_,t_] := (dt)/(dc) * r + (tc)/(dc) * s
In[4]: x[p_,q_,r_,s_,a_,b_,c_,d_,t_] := (ct)/(ca) * u[p,q,r,s,a,b,c,d,t] + (ta)/(ca) * v[p,q,r,s,a,b,c,d,t]
In[5]: y[p_,q_,r_,s_,a_,b_,c_,d_,t_] := (dt)/(db) * v[p,q,r,s,a,b,c,d,t] + (tb)/(db) * w[p,q,r,s,a,b,c,d,t]
In[6]: z[p_,q_,r_,s_,a_,b_,c_,d_,t_] := (ct)/(cb) * x[p,q,r,s,a,b,c,d,t] + (tb)/(cb) * 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
\[a = 0,\] \[b = e,\] \[c = e + f,\] \[d = e+f+g.\]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
\[\frac{d}{dt}Z(P,Q,R,S,0,e,e+f,e+f+g,t)\rvert_{t=e+f} = \frac{f^2 (SR)+g^2 (RQ)}{f g (f+g)}.\]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
\[0,\] \[f,\] \[f+g,\] \[f+g+h.\]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
\[\frac{d}{dt}Z(Q,R,S,T,0,f,f+g,f+g+h,t)\rvert_{t=f} = \frac{f^2 (SR)+g^2 (RQ)}{f g (f+g)}.\]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, CatmullRom 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 nonnegative numbers. Their exact form is irrelevant.
The next step would be to look at the second derivatives, and without further ado we get:
\[\frac{2 \left(f \left(2 e^2 (RS)+3 e g (RQ)+g^2 (PQ)\right)+f^2 (2 e (RS)+g (PQ))e g (2 e+g) (QR)\right)}{e f g (e+f) (f+g)}\]and
\[\frac{2 \left(g \left(f^2 (ST)+3 f h (SR)+2 h^2 (QR)\right)+g^2 (f (ST)+2 h (QR))f h (f+2 h) (RS)\right)}{f g h (f+g) (g+h)}\]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, CatmullRom splines aren’t two times continuously differentiable.
Final note
There’s much more to explore here.
First, the form of the first derivate
\[\frac{f^2 (SR)+g^2 (RQ)}{f g (f+g)}\]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 CatmullRom 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!

Yuksel, C., Schaefer, S., & Keyser, J. (2011). Parameterization and applications of Catmull–Rom curves. ComputerAided Design, 43(7), 747755. ↰

In a better timeline, I would have written a nice overview of CatmullRom splines that I could link to. Alas, here we are. ↰

CITATION NEEDED. ↰

Again, many good reasons I don’t want to get into here. ↰

Apparently, ‘spline’ is the word for curved wooden strips used in building ships and aeroplanes. ↰

Catmull, E., & Rom, R. (1974). A class of local interpolating splines. In Computer aided geometric design (pp. 317326). academic Press. ↰

Barry, P. J., & Goldman, R. N. (1988). A recursive evaluation algorithm for a class of CatmullRom splines. ACM SIGGRAPH Computer Graphics, 22(4), 199204. ↰

The value \(\alpha=\frac{1}{2}\) leads to socalled centripetal CatmullRom splines which are the “best” in some regard.^{1} ↰

‘Lerp’ stands for linear interpolation; the most important function in graphics programming. ↰