Procedural Animation Video Notes

For a first post here, I thought I'd revisit some points that were cut from the previous video, which some viewers asked about in the comments, as well as some common questions.

Isn't all this math expensive to compute?

This is a surprisingly common question, but while the derivations go through a lot of symbols and mathematical machinery, the final implementation is really quite simple, only involving adding and multiplying a handful of numbers per frame.

By comparison, using something like an animation curve, which is a more traditional solution, involves solving a cubic equation every frame, not to mention the significantly larger amount of data that must be stored. Yet even animation curves are cheap enough to be used in all but the most critical cases.

What about FixedUpdate?

FixedUpdate is a good way of equalizing behavior on different devices by locking the update rate to some commonly achievable amount, but by itself, does not guarantee stability. The stability criterion depend on the relationship between the frequency of the system and the frequency of the update, so if we ever want to guarantee stability when the frequency of the system is outside of our control, we need to do the analysis.

Of course, for games that target a high refresh rate, FixedUpdate won't compute the smoothest possible motion, so then you end up needing to do interpolation between the outputs.

Can I apply these dynamics to Quaternions/Orientations?

Quaternions make up a 4-dimensional vector space, so we can use the derivations directly on quaternion objects as long as we implement addition and scalar multiplication. However, the quaternions used for representing rotation are normalized to unit length, and any given orientation has 2 quaternions that describe the same orientation, so there are some pre/postprocessing steps we have to perform:

1. Orient the input quaternion x to face the same direction as state y
2. Compute the state update function as described in the video
3. Normalize the output quaternion y before returning to the user

This isn't a particularly physically-correct way to simulate rotational dynamics, because we are operating on unnormalized quaternions, which do not represent physical rotations. To more correctly simulate rotational dynamics, you would need to use the quaternion operations that preserve normalization. For semi-implicit Euler,

ẋ[n+1] = (x[n+1] - x[n]) / T
y[n+1] = y[n] + Tẏ[n]
ẏ[n+1] = ẏ[n] + T(x[n+1] + k₃ẋ[n+1] - y[n+1] - k₁ẏ[n]) / k₂

becomes something like:

ẋ[n+1] = Log(x[n+1]x⁻¹[n]) / T
y[n+1] = exp(Tẏ[n])y[n]
ẏ[n+1] = ẏ[n] + T(Log(x[n+1]y⁻¹[n+1]) + k₃ẋ[n+1] - k₁ẏ[n]) / k₂

Where Log is the (principal) quaternion logarithm, and exp is the quaternion exponential. However, for small rotations, and in cases where exact physical correctness is not important, I think it's it's fine to use the linear approximations.

Where do your definitions for f, ζ, and r come from?

f, or more commonly, ωₙ=2πf, is the natural frequency of a resonant system, and a typical introduction to linear systems will explore the damped simple harmonic oscillator, which also includes a damping coefficient ζ. This system can be used to model many different types of physical behaviors. The characteristic polynomial of a second-order system is usually presented as the Laplace-domain transfer function:

ωₙ² / (s²+2ζωₙs+ωₙ²)

Which is fine for describing the motion of a mass that's attached to a fixed spring and damper, but if you do the derivation for a mass being dragged by a spring and damper, you instead get the transfer function:

(ωₙζs+ωₙ²) / (s²+2ζωₙs+ωₙ²)

In order to support simulating both of these conditions, I introduced the scaling factor r that scales the location of the zero:

(rωₙζs+ωₙ²) / (s²+2ζωₙs+ωₙ²)

In this way, when r is 0, we get the original equation, and when r is 1, we get the second equation. I realized that we can also place r at other locations, because it affects the initial response of the system, which can be derived using the initial value theorem.

ẏ₀ = lim s→∞ s (rωₙζs+ωₙ²) / (s²+2ζωₙs+ωₙ²)
    = lim s→∞ (rωₙζs²) / (s²+2ζωₙs+ωₙ²) + (ωₙ²s) / (s²+2ζωₙs+ωₙ²)
    = (rωₙζ) / (ωₙ²)
    = rζ / ωₙ

This works nicely because the initial speed, ẏ₀, gets scaled with the natural frequency, meaning changing the natural frequency won't change the shape of the impulse response, so this is the choice of parameterization I decided to go with in the video. It's also nice that the special case of ζ=1, r=1 simplifies to:

   (ωₙs+ωₙ²) / (s²+2ωₙs+ωₙ²)
= ωₙ(s+ωₙ) / (s+ωₙ)²
= ωₙ / (s+ωₙ)

which is a first-order system, so our parameters are also capable of simulating first-order responses (also known as exponential smoothing).

Why Semi-Implicit Euler?

There wasn't any particularly strong reason. I initially went with the velocity Verlet method, but found in analysis that it produces the same discrete-time characteristic equation as semi-implicit Euler, so I switched to the method that was mathematically simpler to explain and computationally simpler to implement.

After the video published, a viewer pointed out that for solving my 2nd order equation, the fully-implicit Euler method produced unconditionally stable results. So rather than:

y[n+1] = y[n] + Tẏ[n]
ẏ[n+1] = ẏ[n] + T(x[n+1] + k₃ẋ[n+1] - y[n+1] - k₁ẏ[n]) / k₂

we have:

y[n+1] = y[n] + Tẏ[n]
ẏ[n+1] = ẏ[n] + T(x[n+1] + k₃ẋ[n+1] - y[n+1] - k₁ẏ[n+1]) / k₂

The second equation has ẏ[n+1] on both sides, so it has to be rearranged in order to compute it:

ẏ[n+1](1 + Tk₁/k₂) = ẏ[n] + T(x[n+1] + k₃ẋ[n+1] - y[n+1]) / k₂
                  ẏ[n+1] = (k₂ẏ[n] + T(x[n+1] + k₃ẋ[n+1] - y[n+1])) / (k₂ + Tk₁)

And if we perform the stability analysis on this discrete-time system, we can see the eigenvalues never get bigger than a magnitude of 1.

This is what I use in my game now, but I still keep the pole matching solution when ωₙT > ζ because it produces more accurate results when taking large steps, or simulating very stiff systems. That reminds me...

What's Pole Matching?

Pole matching is a method of discretizing a differential equation for computation that's slightly better than numerical integration techniques, but isn't too expensive to re-compute in real-time like some other methods. Instead of using the difference equation from the pole-matching technique directly, I use it to place the z-domain poles of the Euler method discretization by computing correction terms for k₁ and k₂.

Essentially, I use pole-matching to find the theoretical location of the z-domain poles given T, k₁, k₂, and then derive how much I would need to modify k₁ and k₂ such that the z-domain poles of the Euler method system get put there. if the sampling period T is variable, or if changes are made to k₁ or k₂, then this modification will need to be recomputed as needed.

The reason for this whole dance is so that the state variables (y and ẏ) remain unchanged, allowing for seamless switching between regular operation and pole-matching operation.