Perfect Tracking with Springs

Created on Jan. 25, 2023, 3:44 p.m.

If you have some animation data with discontinuities you want to remove, you might have considered tracking that animation data with some kind of spring:

This does a good job at removing the discontinuities - but the spring needs to be very stiff - meaning the discontinuities get aggressively blended out - and even worse - it never quite reproduces exactly the input signal in places where nothing is wrong.

For this reason, I've always tried to avoid going down this route whenever I could, and to me, inertialization always seemed like the "correct" way of dealing with these kind of discontinuities if you knew where they occurred.

However, a few years ago, my old colleague Jack Potter showed me a kind of spring that can filter animation in exactly the way we want: it can remove discontinuities smoothly while still perfectly tracking the original signal when nothing is wrong.

The trick is to make a spring that tracks the acceleration, velocity, and position of the input animation in different proportions.

In code it looks something like this:

void tracking_spring_update(
    float& x,
    float& v,
    float x_goal,
    float v_goal,
    float a_goal,
    float x_gain,
    float v_gain,
    float a_gain,
    float dt)
{
    v = lerp(v, v + a_goal * dt, a_gain);
    v = lerp(v, v_goal, v_gain);
    v = lerp(v, (x_goal - x) / dt, x_gain);
    x = x + dt * v;
}

First, we blend the current velocity with the current velocity plus the acceleration of the input animation multiplied by the dt. Then, we blend this with the velocity of the input animation itself. Finally, we blend this with a velocity that pulls us toward the position of the input animation.

For each of these blends we can use a gain to control the strength. For perfect tracking, we will want to set the acceleration gain to 1, the velocity gain to something like 0.2, and the position gain to something quite small such as 0.01 (when running at 60 frames per second).

This spring will need to be fed the acceleration, velocity, and position of the input signal. The easiest way to get these is by finite difference:

float tracking_target_acceleration(
    float x_next,
    float x_curr,
    float x_prev,
    float dt)
{
    return (((x_next - x_curr) / dt) - ((x_curr - x_prev) / dt)) / dt;
}

float tracking_target_velocity(
    float x_next,
    float x_curr,
    float dt)
{
    return (x_next - x_curr) / dt;
}

If we know there is a discontinuity in this signal, we know the values computed by this finite difference will not be good (at least the acceleration and velocity targets wont be) - so we just ignore them for those time steps and only blend the velocity with values we know are good:

void tracking_spring_update_no_acceleration(
    float& x,
    float& v,
    float x_goal,
    float v_goal,
    float x_gain,
    float v_gain,
    float dt)
{
    v = lerp(v, v_goal, v_gain);
    v = lerp(v, (x_goal - x) / dt, x_gain);
    x = x + dt * v;
}

void tracking_spring_update_no_velocity_acceleration(
    float& x,
    float& v,
    float x_goal,
    float x_gain,
    float dt)
{
    v = lerp(v, (x_goal - x) / dt, x_gain);
    x = x + dt * v;
}

This allows the spring to ignore the discontinuity and naturally converge back onto the input stream, tracking it perfectly:

If we don't know where the discontinuities are in the input signal we have two options - we can try to detect them and ignore them based on some kind of heuristic - or we can clip whatever we get to some maximum and minimum velocity and acceleration. Clipping can still sometimes give us a nasty jump when discontinuities occur but works okay given the imperfect situation we might be in:

This can be useful when our input signal may jump around in hard to specify ways or when it is coming from some black box we can't control.

One important thing to note here is that this code is not in any way robust to varying timesteps! The use of these gains will give us very different results if we tick at different rates:

We can improve things somewhat by switching lerp to damper_exact and using half-lives to control the blending (in this case I use 0.0 for the acceleration halflife, 0.05 for the velocity halflife and 1.0 for the position halflife).

void tracking_spring_update_improved(
    float& x,
    float& v,
    float x_goal,
    float v_goal,
    float a_goal,
    float x_halflife,
    float v_halflife,
    float a_halflife,
    float dt)
{
    v = damper_exact(v, v + a_goal * dt, a_halflife, dt);
    v = damper_exact(v, v_goal, v_halflife, dt);
    v = damper_exact(v, (x_goal - x) / dt, x_halflife, dt);
    x = x + dt * v;
}

This version is certainly far from dt-invariant either, but feels a bit better than before to me.

Although basically stable, I would warn against ticking this spring at a variable timestep. Other than that I think it is pretty neat.

If anyone wants to go through the maths to try and derive an exact, dt-invariant version of this spring that would be awesome - or if anyone knows if this spring exists under a different name in engineering/control theory I would love to know!

That's all for now.

Appendix: Exact Tracking Spring

So after a little discussion on Twitter it seemed that (although there were some good suggestions), no one could identify exactly what kind of spring this was. A couple of the replies however did seed some ideas in my head and I decided to try and see if I could derive a dt-invariant version of this spring by playing around with the equations. Spoiler: it turns out this spring is equivalent to our normal spring-damper with specifically tuned stiffness, damping, goal position, and goal velocity.

To see how this works, we need to break down the tracking spring velocity equation into its different parts.

First, we have the continuation of the existing velocity, which will end up scaled by some amount depending on our gains. We will denote this scaling using some variable \( \theta_0 \):

\begin{align*} v_{t} &= \theta_0 \cdot v_t + ... \end{align*}

Next, we have the contribution from the acceleration of the signal being tracked \( g_a \). This is multiplied by the delta time \( dt \), and again scaled by some amount depending on the gains. Similarly to before, we'll denote this as \( \theta_1 \):

\begin{align*} v_{t} &= \theta_0 \cdot v_t + dt \cdot \theta_1 \cdot g_a + ... \end{align*}

Then, we have the contribution from the velocity of the signal being tracked \( g_v \), which we will multiply by the delta time \( dt \) and our third scaling variable \( \theta_2 \):

\begin{align*} v_{t} &= \theta_0 \cdot v_t + dt \cdot \theta_1 \cdot g_a + dt \cdot \theta_2 \cdot g_v + ... \end{align*}

Finally, we have the term that pulls the spring toward the position of the signal being tracked \( g_x \). Like the others, this we will scale by the delta time \( dt \) and another variable \( \theta_3 \):

\begin{align*} v_{t} &= \theta_0 \cdot v_t + dt \cdot \theta_1 \cdot g_a + dt \cdot \theta_2 \cdot g_v + dt \cdot \theta_3 \cdot (g_x - x_{t}) \end{align*}

We can think of this as simply another formulation of our tracking spring using a different set of variables. What we need to do now is find a way to compute these values \( \theta_0 \), \( \theta_1 \), \( \theta_2 \), \( \theta_3 \) from our gains. The way we are going to do this is to kind of work backwards to find assignments of these that reflect the behavior we expect from the lerp formulation.

First, what are the values of \( \theta_0 \) and \( \theta_1 \) according to our acceleration gain \( \alpha_a \)? This one is simple - the acceleration gain \( \alpha_a \) doesn't affect the scale of the existing velocity, meanwhile when \( \alpha_a = 1 \) we expect to have the full contribution from \( g_a \) - effectively \( \theta_1 = \alpha_a \) in this case:

\begin{align*} \theta_0 &= 1 \\ \theta_1 &= \alpha_a \\ \end{align*}

Next, what happens when we introduce the velocity gain \( \alpha_v \)? As well as scaling our velocity target \( g_v \) up such that \( dt \cdot \theta_2 = 1 \) when \( \alpha_v = 1 \), it will also scale down the values of \( \theta_0 \) and \( \theta_1 \) in equal proportion.

\begin{align*} \theta_0 &= 1 - \alpha_v \\ \theta_1 &= \alpha_a \cdot (1 - \alpha_v) \\ \theta_2 &= \frac{\alpha_v}{dt} \end{align*}

Finally, we need to account for our position gain \( \alpha_x \). When \( \alpha_x = 1 \) we expect the result to be pulled directly onto the tracking target \( g_x \), which means \( \alpha_3 = \tfrac{\alpha_x}{dt^2} \). Also, just like the velocity gain, the position gain will scale down the other terms as it grows larger:

\begin{align*} \theta_0 &= (1 - \alpha_v) \cdot (1 - \alpha_x) \\ \theta_1 &= \alpha_a \cdot (1 - \alpha_v) \cdot (1 - \alpha_x) \\ \theta_2 &= \frac{\alpha_v \cdot (1 - \alpha_x)}{dt} \\ \theta_3 &= \frac{\alpha_x}{dt^2} \end{align*}

But why is that useful? Well our goal now is going to be to rearrange our equation with \( \theta_0 \), \( \theta_1 \), \( \theta_2 \) and \( \theta_3 \) to put it in a similar form to our spring-damper velocity equation. If we can find a matching formulation then we can use the time-invariant derivation we've already worked out.

First let's re-arrange our spring-damper velocity equation a bit, where \( s = stiffness \), \( d = damping \), \( g \) is our goal position, and \( q \) is our goal velocity:

\begin{align*} v_t = v_t + dt \cdot s \cdot (g - x_t) + dt \cdot d \cdot (q - v_t) \end{align*}

First we'll multiply out the right-most term:

\begin{align*} v_t = v_t + dt \cdot s \cdot (g - x_t) + dt \cdot d \cdot q - dt \cdot d \cdot v_t \end{align*}

Then, we can re-factorize to put the single velocity term along with the rightmost term to get this:

\begin{align*} v_t = dt \cdot s \cdot (g - x_t) + dt \cdot d \cdot q + (1 - dt \cdot d) \cdot v_t \end{align*}

Now, let's take our other equation...

\begin{align*} v_{t} &= \theta_0 \cdot v_t + dt \cdot \theta_1 \cdot g_a + dt \cdot \theta_2 \cdot g_v + dt \cdot \theta_3 \cdot (g_x - x_{t}) \end{align*}

... and factorize out the terms containing \( \theta_1 \) and \( \theta_2 \) so they share the same \( dt \):

\begin{align*} v_{t} &= \theta_0 \cdot v_t + dt \cdot \theta_1 \cdot g_a + dt \cdot \theta_2 \cdot g_v + dt \cdot \theta_3 \cdot (g_x - x_{t}) \\ v_{t} &= dt \cdot \theta_3 \cdot (g_x - x_{t}) + dt \cdot (\theta_1 \cdot g_a + \theta_2 \cdot g_v) + \theta_0 \cdot v_t \end{align*}

Aha! If we look carefully we can see some matches. For the first term we have \( s = \theta_3 \), and \( g = g_x \). For the second term we have \( d \cdot q = \theta_1 \cdot g_a + \theta_2 \cdot g_v \). And for the third term we have \( 1 - dt \cdot d = \theta_0 \). While the first two equations give us direct solutions, the second two will require a little re-arranging as we need to solve for \( d \) and \( q \).

We can solve for \( d \) first:

\begin{align*} 1 - dt \cdot d &= \theta_0 \\ dt \cdot d &= 1 - \theta_0 \\ d &= \frac{1 - \theta_0}{dt} \end{align*}

Which we can then use to solve for \( q \)

\begin{align*} d \cdot q &= \theta_1 \cdot g_a + \theta_2 \cdot g_v \\ q &= \frac{\theta_1 \cdot g_a + \theta_2 \cdot g_v}{d} \end{align*}

And with that we have our spring-damper parameters \( s \), \( d \), \( g \), \( q \), computed from our tracking spring parameters \( \theta_0 \), \( \theta_1 \), \( \theta_2 \) and \( \theta_3 \):

\begin{align*} s &= \theta_3 \\ d &= \frac{1 - \theta_0}{dt} \\ g &= g_x \\ q &= \frac{\theta_1 \cdot g_a + \theta_2 \cdot g_v}{d} \\ \end{align*}

Now all that remains is to plug these parameters into the exact spring damper equation we already derived and we are done. Here it is implemented in C++, where gain_dt is the timestep which we are emulating running our gains at:

void tracking_spring_update_exact(
    float& x,
    float& v,
    float x_goal,
    float v_goal,
    float a_goal,
    float x_gain,
    float v_gain,
    float a_gain,
    float dt,
    float gain_dt)
{
    float t0 = (1.0f - v_gain) * (1.0f - x_gain);
    float t1 = a_gain * (1.0f - v_gain) * (1.0f - x_gain);
    float t2 = (v_gain * (1.0f - x_gain)) / gain_dt;
    float t3 = x_gain / (gain_dt*gain_dt);
    
    float stiffness = t3;
    float damping = (1.0f - t0) / gain_dt;
    float spring_x_goal = x_goal;
    float spring_v_goal = (t2*v_goal + t1*a_goal) / ((1.0f - t0) / gain_dt);
    
    spring_damper_exact_stiffness_damping(
      x, 
      v, 
      spring_x_goal,
      spring_v_goal,
      stiffness,
      damping,
      dt);
}

And this is what it looks like in action.

In terms of dt-invariance definitely much better than what we had before!

However, one thing you've probably noticed about this version is that it no longer tracks the input signal pixel perfectly. This is because the call to spring_damper_exact does not perfectly emulate what you would get evaluating the spring-damper with discrete ticks at a rate of gain_dt - which is what we are assuming when we compute our spring parameters. I think if we were to set gain_dt to a smaller value and tweak the gains we would get better tracking. Alternatively, we could switch to the normal integration based method when running at the framerate specified by gain_dt.

Another thing to note is that this derivation assumes that the target position, velocity and acceleration coming from the input signal are fixed over the full timestep. This means that if you make the \( dt \) very large and your signal becomes very coarse the behavior of this spring will be different to what you would get with smaller timesteps. In other words, this dt-invariant derivation becomes less and less effective if the tracking signal is changing quickly.

Nonetheless, if you do have a timestep that varies, I think this version is still is far, far better than relying on the version we had before.

The other variations we need where we don't have acceleration or velocity targets have pretty similar derivations so I wont go through the maths, but in C++ they can be implemented as follows:

void tracking_spring_update_no_acceleration_exact(
    float& x,
    float& v,
    float x_goal,
    float v_goal,
    float x_gain,
    float v_gain,
    float dt,
    float gain_dt)
{
    float t0 = (1.0f - v_gain) * (1.0f - x_gain);
    float t2 = (v_gain * (1.0f - x_gain)) / gain_dt;
    float t3 = x_gain / (gain_dt*gain_dt);
    
    float stiffness = t3;
    float damping = (1.0f - t0) / gain_dt;
    float spring_x_goal = x_goal;
    float spring_v_goal = t2*v_goal / ((1.0f - t0) / gain_dt);

    spring_damper_exact_stiffness_damping(
      x, 
      v, 
      spring_x_goal,
      spring_v_goal,
      stiffness,
      damping,
      dt);
}

void tracking_spring_update_no_velocity_acceleration_exact(
    float& x,
    float& v,
    float x_goal,
    float x_gain,
    float dt,
    float gain_dt)
{
    float t0 = 1.0f - x_gain;
    float t3 = x_gain / (gain_dt*gain_dt);
    
    float stiffness = t3;
    float damping = (1.0f - t0) / gain_dt;
    float spring_x_goal = x_goal;
    float spring_v_goal = 0.0f;
  
    spring_damper_exact_stiffness_damping(
      x, 
      v, 
      spring_x_goal,
      spring_v_goal,
      stiffness,
      damping,
      dt);
}

And that is about it! All the code for this article can be found on github. I hope you find it useful.