# A Cubic Spline for Animation

I am designing animations for a string of addressable LEDs. Deterministic patterns generally have a nice regularity to them, but their visual interest diminishes after a while. Random patterns can have a bit more visual interest, since they are different every time. If we just assign a random value to each pixel, the result is very busy, loud, disjointed. It would be better to have something changing in a random way, but with a smooth, regular structure so that the eye can actually take it all in.

One nice way to get a smooth, organic-feeling function is to use a polynomial spline. We can specify values to achieve at certain times (in animation these are called "key frames") and use a polynomial to interpolate between them. By choosing the spline carefully, we can get a motion which is not only continuous, but smooth in both time and space. It will have no jumps, and no times or places where the color suddenly starts changing in a different way. In this post we will walk through a procedure for generating those splines using linear algebra.

# The setting

We have a one-dimensional string of LEDs. We can specify the location of each LED by a number on an interval. To make the numbers work out nice, let's put them on the interval $$[0, 3]$$. For symmetry, let's describe time in the same way, as an interval $$[0, 3]$$. For reasons which will become clear, we will focus on the time subinterval $$[1, 2]$$. We will put the key frames at $$t=0$$, $$t=1$$, $$t=2$$, and $$t=3$$. We will specify each key frame by the values at $$s=0$$, $$s=1$$, $$s=2$$, and $$s=3$$. That is, we begin with a function defined on $$\{0, 1, 2, 3\}\times\{0, 1, 2, 3\}$$, and we wish to extend it to a function on $$[1, 2]\times[0, 3]$$.

Why use this specific number of key values? I want to make a cubic polynomial, which we can write as $$f(t, s) = \sum_{i, j=0}^3 A_{ij}t^i s^j$$. This has sixteen degrees of freedom. By specifying the function at sixteen points, we can make a system which is neither over-constrained nor under-constrained.

Why go for a cubic function? I want to get a curve which is smooth in time and space. On a given interval, I want it to hit the key frames, but I also want it to be smooth at both ends. That makes four constraints, so a third degree polynomial is just right.

Why look at times between 1 and 2? I want to think of time as starting at 1, and proceeding to 2. When the time hits 2, we can extend $$f$$ by defining $$f(t, s) = g(t-1, s)$$ for $$t\in\{1\}\cup[2,3]$$. In this way, $$g$$ has the same form as $$f$$, with values of $$g(t)$$ specified for $$t\in \{0, 1, 2\}$$. By dropping the values of $$f$$ at $$t=0$$ and making new values for $$g$$ at $$t=3$$, the program repeats with $$g$$ taking the role of $$f$$.

# The clever bit

When we restrict $$f$$ to $$t=2$$, we get a polynomial in $$s$$. Since it interpolates four values at $$s=0, 1, 2, 3$$, it is completely specified. Since $$g$$ restricted to $$t=1$$ is also a cubic interpolant, we can see that $$f$$ is continuous not only at $$\{2\}\times \{0, 1, 2, 3\}$$ but on all of $$\{2\}\times [0, 3]$$. If we are clever about how we specify $$f$$, we can do even better. Notice that $$f_t$$ is also a cubic polynomial in $$s$$ when restricted to $$t=2$$. If we can specify $$f_t(2, s)=g_t(1, s)$$ then we can make $$f'$$ continuous as well!

How to do this? Since $$f$$ is never directly evaluated at $$t=3$$, we will use the funciton values there to instead specify $$f_t$$ at $$t=2$$. That is, we will enforce $$f_t(2, s) = (1-c) (f(3, s)-f(1, s))/2$$, where $$c$$ is a "tension" term allowing us to tweak the behavior of the interpolant. By assigning the same rule for $$g$$, that is, $$g_t(1, s) = (1-c)(f(3, s)-f(1, s))/2$$, we get a function which is continuous and smooth.

# Computing it

We began by writing $$f(t, s)=\sum_{i, j=0}^3 A_{ij}t^i s^j$$. We can write this as a matrix product,

$$f(t, s) = \begin{bmatrix} 1 \\ t \\ t^2 \\ t^3 \end{bmatrix}^T A \begin{bmatrix} 1 \\ s \\ s^2 \\ s^3 \end{bmatrix}.$$

to keep things tidy, let's define $$P_x = [1, x, x^2, x^3]^T$$. Then we can write this as $$f(t, s)=P_t^T A P_s$$. The goal then is to specify the matrix $$A$$, in terms of the data, the key-frame values of $$f$$ on $$\{0, 1, 2, 3\}^2$$. Write $$F=[f(i, j)]_{i, j=0}^3$$, so we are looking for how to write $$A$$ in terms of $$F$$.

## Continuity

We use the values at $$t=1, 2$$ to enforce the continuity of the spline. This is straightforward: we want $$P_1^TAP_s = F_{1s}$$ for $$s=0, 1, 2, 3$$. Writing $$P_S = [P_0P_1P_2P_3]$$ we can write this more compactly as $$P_1^TAP_S = F_1$$. Similarly, we want $$P_2^TAP_S = F_2$$. Together, these give

$$[P_1, P_2]^T A P_S = \begin{bmatrix} 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \end{bmatrix} F.$$

Not enough to solve for $$A$$ yet, but we have not used smoothness yet.

## Smoothness

Differentiation yields

$$f_t(t, s) = \begin{bmatrix} 0 \\ 1 \\ 2t \\ 3t^2 \end{bmatrix}^T A \begin{bmatrix} 1 \\ s \\ s^2 \\ s^3 \end{bmatrix}.$$

Write $$D_x=\frac{d}{dx} P_x$$, so we can write this as $$f_t = D_t^T A P_s$$. We want $$f(1, s) = \frac{1-c}{2}(f(2, s) - f(0, s))$$. As a matrix equation, this becomes

$$D_1^T A P_S = \frac{1-c}{2} \begin{bmatrix}-1 & 0 & 1 & 0\end{bmatrix}F.$$

Repeating for the $$t=2$$ case yields

$$[D_1, D_2]^T A P_S = \frac{1-c}{2}\begin{bmatrix}-1&0&1&0\\0&-1&0&1\end{bmatrix}F.$$

Writing it this way shows how to combine with the continuity equation,

$$[P_1, P_2, D_1, D_2]^T A P_S = \frac{1}{2}\left( \begin{bmatrix} 0 & 2 & 0 & 0\\ 0 & 0 & 2 & 0\\ -1& 0 & 1 & 0\\ 0 &-1 & 0 & 1 \end{bmatrix} +c\begin{bmatrix} 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ 1 & 0 &-1 & 0\\ 0 & 1 & 0 &-1 \end{bmatrix} \right)F.$$

## Solving for $$A$$

Now all the matrices are square, so we are in business. We have enough to solve for $$A$$. Writing $$U=[P_1,P_2,D_1,D_2]^T$$ we find

$$U^{-1} = \begin{bmatrix}-4 & 5 & -4 & -2\\12 & -12 & 8 & 5\\-9 & 9 & -5 & -4\\2 & -2 & 1 & 1\end{bmatrix}.$$

We can also compute

$$P_S^{-1} = \begin{bmatrix}1 & 1 & 1 & 1\\0 & 1 & 2 & 3\\0 & 1 & 4 & 9\\0 & 1 & 8 & 27\end{bmatrix}^{-1} =\frac{1}{6} \begin{bmatrix}6 & -11 & 6 & -1\\0 & 18 & -15 & 3\\0 & -9 & 12 & -3\\0 & 2 & -3 & 1\end{bmatrix}.$$

Multiplying on both sides yields

$$A = \frac{1}{12} \left(\begin{bmatrix}4 & -6 & 6 & -2\\-8 & 19 & -16 & 5\\5 & -14 & 13 & -4\\-1 & 3 & -3 & 1\end{bmatrix} + c \begin{bmatrix}-4 & -2 & 4 & 2\\8 & 5 & -8 & -5\\-5 & -4 & 5 & 4\\1 & 1 & -1 & -1\end{bmatrix} \right) F \begin{bmatrix}-4 & 5 & -4 & -2\\12 & -12 & 8 & 5\\-9 & 9 & -5 & -4\\2 & -2 & 1 & 1\end{bmatrix}.$$

## Applying the formula

Now that we have a formula for $$f(t, s)=P_t^TAP_s$$, we will want to apply it to our light string, which is discrete. Choose some discrete times $$T = [1, 1+\Delta t, 1+2\Delta t, \dots, 2]$$ and string positions $$S=[0,\Delta s, 2\Delta s, ..., 3]$$. Form them into polynomials $$P_T$$ and $$P_S$$ as before, then multiply them on the outside in advance; the coefficient matrices we derived above will not change, only F will change. This can save us a lot of multiplication, since much of it can be done in advance.

# Conclusion

I won't belabor this any further here. You can see how this algorithm is implemented in the repository I have shared. I think this shows how linear algebra allows us to express some pretty complicated formulas quite compactly, which lets us think about more-complicated things than would otherwise be possible.