I'm trying to reverse engineer an interpolation used for animation. It seems to use 5 coefficients of the taylor series of $f(x) =e^{x}$.
import math
import numpy as np
def exp_taylor_series(x):
return [x**n / math.factorial(n) for n in range(5)]
def interp(t, p0, p1, v0, v1, t0, t1):
"""
t : The time at which to calculate the interpolated value.
p0 : The position of the first point.
p1 : The position of the second point.
v0 : The velocity at the first point.
v1 : The velocity at the second point.
t0 : The time at the first point.
t1 : The time at the second point.
Returns: The interpolated value of the point at time t.
"""
# Get x at known point p1 with slope m1 at time t1
# Use known point p0 with slope m0 at time t0 as offset
p = p1 - p0 - v0 * (t1 - t0)
v = v1 - v0
y = np.array([p, v, 0.0])
s = exp_taylor_series(t1 - t0)
A = np.array([
[s[4], s[3], s[2]],
[s[3], s[2], s[1]],
[s[2], s[1], s[0]],
])
x = np.dot(np.linalg.inv(A), y) # x = A^-1 ⋅ y
# Do the inverse of previous
# Solve for y at current time t to get slope and position
# Use known point p0 with slope m0 at time t0 as offset
s = exp_taylor_series(t - t0)
A = np.array([
[s[4], s[3], s[2]],
[s[3], s[2], s[1]],
[s[2], s[1], s[0]],
])
y = np.dot(A, x) # y = A ⋅ x
p = y[0]
v = y[1]
pt = p + p0 + v0 * (t - t0)
vt = v + v0
return pt
For $p_0 = 0, p_1 = 1, m_0 = 0, m_1 = 0, t_0 = 0, t_1 = 1$, it gives a curve that looks a lot like a Hermite Spline.
Something notable that I saw is that the difference between the two curves is a quadradic equation with the form of $d(x) = 2 - 4x + 2x^2 + 5.95499956\times10^{-15}x -3.53508224\times10^{-16}x $
What is this matrix supposed to represent and does it have a name?
Edit
The curve is in fact a scaled and translated quintic curve which passes through $(0, 0)$ with a slope of $m_0$ and through $(1, 1)$ with a slope of $m_1$.
I did some testing because I wondered if normalizing the curve could reduce the number of variables at play. To achieve this, I made a different function where $t \in [0, 1]$ and $p_0 = 0, p_1=1$ leaving only the start and end velocities as free variables and they can be seen as slopes now.
I was able to reduce the function to a simple matrix multiplication and get the same curve.
The matrix comes from starting with the taylor series expansion of exp(x) at 1 $$ s(x) = \left(\frac{x^0}{0!}, \frac{x^1}{1!}, ..., \frac{x^n}{n!}, ...\right) $$ Constructing the a matrix to satisfy $y = Ax$ using the first 5 terms for solving a linear equation. The construction is odd to me. I do not know why it's done this way. $$ A = \begin{bmatrix} s_4 & s_3 & s_2\\ s_3 & s_2 & s_1\\ s_2 & s_1 & s_0 \end{bmatrix} = \begin{bmatrix} \frac{1}{24} & \frac{1}{6} & \frac{1}{2} \\ \frac{1}{6} & \frac{1}{2} & 1 \\ \frac{1}{2} & 1 & 1 \end{bmatrix} $$
Since this matrix construction is constant and no longer dependent on different values of $t_0$ and $t_1$, the expensive computation of the inverse can be calculated and stored.
$$ A^{-1} = \begin{bmatrix} 72 & -48 & 12 \\ -48 & 30 & -6 \\ 12 & -6 & 1 \end{bmatrix} $$
We normally get y by using range values for position, velocity and acceleration. However, thanks to our assumptions this is greatly simplified.
$$ y = \begin{bmatrix} p_1 - p_0 - m_0 * (t_1 - t_0) \\ m_1 - m_0 \\ 0 \end{bmatrix} = \begin{bmatrix} 1 - m_0 & m_1 - m_0 & 0 \end{bmatrix} $$
Now our formula for $x$ is
$$ x = \begin{bmatrix} 72 & -48 & 12 \\ -48 & 30 & -6 \\ 12 & -6 & 1 \end{bmatrix} \begin{bmatrix} 1 - m_0 & m_1 - m_0 & 0 \end{bmatrix} $$
At this point the expensive bit of computing the inverse matrix for the linear equation is done and the interpolation can now use it to get new values to x with $y = A_t x$. Where $A_t$ is constructed again with $s(x)$ at $t$.
$$ y_t = A_t x = \begin{bmatrix} \frac{t^4}{24} & \frac{t^3}{6} & \frac{t^2}{2} \\ \frac{t^3}{6} & \frac{t^2}{2} & t \\ \frac{t^2}{2} & t & 1 \end{bmatrix} \begin{bmatrix} 72 & -48 & 12 \\ -48 & 30 & -6 \\ 12 & -6 & 1 \end{bmatrix} \begin{bmatrix} 1 - m_0 & m_1 - m_0 & 0 \end{bmatrix} $$
I then reconstruct the position, velocity and acceleration from $y_t$ using the reverse process. $$ p_t = y_{t_0} + p_0 + m_0 \dot (t - t_0) = y_{t_0} + m_0 t $$ $$ m_t = y_{t_1} + m_0 $$
I don't have the work for this, but using the polyval function in python, I can create a degree 5 polynomial curve and this lets me skip finding $y_t$ and directly get the first three coefficients to which I add $m_0$ and $p_0 = 0$.
$$ c = \begin{bmatrix} \frac{1^4}{4!} & 0 & 0 \\ 0 & \frac{1^3}{3!} & 0 \\ 0 & 0 & \frac{1^2}{2!} \end{bmatrix} x = \begin{bmatrix} \frac{1}{24} & 0 & 0 \\ 0 & \frac{1}{6} & 0 \\ 0 & 0 & \frac{1}{2} \end{bmatrix} \begin{bmatrix} 72 & -48 & 12 \\ -48 & 30 & -6 \\ 12 & -6 & 1 \end{bmatrix} \begin{bmatrix} 1 - m_0 & m_1 - m_0 & 0 \end{bmatrix} = \begin{bmatrix} 3 & -8 & 6 \\ -2 & 5 & -3 \\ 0.5 & -1 & 0.5 \end{bmatrix} \begin{bmatrix} 1 - m_0 & m_1 - m_0 & 0 \end{bmatrix} $$
Instead of adding the coefficients after the fact, I re-arange the matrix and x to look like:
$$ c = \begin{bmatrix} -2 & -1 & 3 \\ 5 & 3 & -8 \\ -3 & -3 & 6 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \\ \end{bmatrix} \begin{bmatrix} m_1 & m_0 & 1 \end{bmatrix} $$
Which is enough to model the curve. The simplified code looks like this
import numpy as np
def interp_simple(m0, m1, t):
"""
Interpolate for a Quintic curve which passes through (0, 0) with a slope of m0
and through (1, 1) with a slope of m1.
Acceleration, jerk and snap are hardcoded and come from a strange taylor series of the exp function
t : Interpolation value between 0 and 1.
m0 : The slope at t=0.
m1 : The slope at t=1.
Returns: The interpolated value for value t.
"""
M = np.array([
[-2., -1., 3.],
[ 5., 3., -8.],
[-3., -3., 6.],
[ 0., 1., 0.],
[ 0., 0., 0.]
])
return np.polyval(np.dot(M, np.array([m1, m0, 1.0])), t)
def interp_simplified(t, p0, p1, v0, v1, t0, t1):
"""
t : The time at which to calculate the interpolated value.
p0 : The position of the first point.
p1 : The position of the second point.
v0 : The velocity at the first point.
v1 : The velocity at the second point.
t0 : The time at the first point.
t1 : The time at the second point
Returns: The interpolated value of the point at time t.
"""
# Catch a few singularities
if p0 == p1:
return p0
if t0 == t1:
return (p0 + p1) / 2
# Get the ranges
pr = p1 - p0
tr = t1 - t0
vr = pr / tr
# Interpolate like if p0=0, p1=1, t0=0 and t1=1, duration will be 0
# Also determine the slopes from the speeds
return interp_simple(v0 / vr, v1 / vr, t / tr) * pr + p0
The curve fits perfectly with these coefficients
I still would like to know why the taylor series for $e^x$ is used as it would let me know how the higher order coefficients are used.