Icon BVHView


Icon Dead Blending Node in Unreal Engine


Icon Propagating Velocities through Animation Systems


Icon Cubic Interpolation of Quaternions


Icon Dead Blending


Icon Perfect Tracking with Springs


Icon Creating Looping Animations from Motion Capture


Icon My Favourite Things


Icon Inertialization Transition Cost


Icon Scalar Velocity


Icon Tags, Ranges and Masks


Icon Fitting Code Driven Displacement


Icon atoi and Trillions of Whales


Icon SuperTrack: Motion Tracking for Physically Simulated Characters using Supervised Learning


Icon Joint Limits


Icon Code vs Data Driven Displacement


Icon Exponential Map, Angle Axis, and Angular Velocity


Icon Encoding Events for Neural Networks


Icon Visualizing Rotation Spaces


Icon Spring-It-On: The Game Developer's Spring-Roll-Call


Icon Interviewing Advice from the Other Side of the Table


Icon Saguaro


Icon Learned Motion Matching


Icon Why Can't I Reproduce Their Results?


Icon Latinendian vs Arabendian


Icon Machine Learning, Kolmogorov Complexity, and Squishy Bunnies


Icon Subspace Neural Physics: Fast Data-Driven Interactive Simulation


Icon Software for Rent


Icon Naraleian Caterpillars


Icon The Scientific Method is a Virus


Icon Local Minima, Saddle Points, and Plateaus


Icon Robust Solving of Optical Motion Capture Data by Denoising


Icon Simple Concurrency in Python


Icon The Software Thief


Icon ASCII : A Love Letter


Icon My Neural Network isn't working! What should I do?


Icon Phase-Functioned Neural Networks for Character Control


Icon 17 Line Markov Chain


Icon 14 Character Random Number Generator


Icon Simple Two Joint IK


Icon Generating Icons with Pixel Sorting


Icon Neural Network Ambient Occlusion


Icon Three Short Stories about the East Coast Main Line


Icon The New Alphabet


Icon "The Color Munifni Exists"


Icon A Deep Learning Framework For Character Motion Synthesis and Editing


Icon The Halting Problem and The Moral Arbitrator


Icon The Witness


Icon Four Seasons Crisp Omelette


Icon At the Bottom of the Elevator


Icon Tracing Functions in Python


Icon Still Things and Moving Things


Icon water.cpp


Icon Making Poetry in Piet


Icon Learning Motion Manifolds with Convolutional Autoencoders


Icon Learning an Inverse Rig Mapping for Character Animation


Icon Infinity Doesn't Exist


Icon Polyconf


Icon Raleigh


Icon The Skagerrak


Icon Printing a Stack Trace with MinGW


Icon The Border Pines


Icon You could have invented Parser Combinators


Icon Ready for the Fight


Icon Earthbound


Icon Turing Drawings


Icon Lost Child Announcement


Icon Shelter


Icon Data Science, how hard can it be?


Icon Denki Furo


Icon In Defence of the Unitype


Icon Maya Velocity Node


Icon Sandy Denny


Icon What type of Machine is the C Preprocessor?


Icon Which AI is more human?


Icon Gone Home


Icon Thoughts on Japan


Icon Can Computers Think?


Icon Counting Sheep & Infinity


Icon How Nature Builds Computers


Icon Painkillers


Icon Correct Box Sphere Intersection


Icon Avoiding Shader Conditionals


Icon Writing Portable OpenGL


Icon The Only Cable Car in Ireland


Icon Is the C Preprocessor Turing Complete?


Icon The aesthetics of code


Icon Issues with SDL on iOS and Android


Icon How I learned to stop worrying and love statistics


Icon PyMark


Icon AutoC Tools


Icon Scripting xNormal with Python


Icon Six Myths About Ray Tracing


Icon The Web Giants Will Fall


Icon PyAutoC


Icon The Pirate Song


Icon Dear Esther


Icon Unsharp Anti Aliasing


Icon The First Boy


Icon Parallel programming isn't hard, optimisation is.


Icon Skyrim


Icon Recognizing a language is solving a problem


Icon Could an animal learn to program?




Icon Pure Depth SSAO


Icon Synchronized in Python


Icon 3d Printing


Icon Real Time Graphics is Virtual Reality


Icon Painting Style Renderer


Icon A very hard problem


Icon Indie Development vs Modding


Icon Corange


Icon 3ds Max PLY Exporter


Icon A Case for the Technical Artist


Icon Enums


Icon Scorpions have won evolution


Icon Dirt and Ashes


Icon Lazy Python


Icon Subdivision Modelling


Icon The Owl


Icon Mouse Traps


Icon Updated Art Reel


Icon Tech Reel


Icon Graphics Aren't the Enemy


Icon On Being A Games Artist


Icon The Bluebird


Icon Everything2


Icon Duck Engine


Icon Boarding Preview


Icon Sailing Preview


Icon Exodus Village Flyover


Icon Art Reel




Icon One Cat Just Leads To Another

Inertialization Transition Cost

Created on June 2, 2022, 3:21 p.m.

Something which is often required by animation systems (Motion Matching being a key example) is a way to compute a "cost" associated with a particular transition between two frames of animation.

In Motion Matching this is typically done by taking the difference between "features" of both the source and destination animations. Usually, the positions of a few key joints (such as the feet), as well as their velocities. The magnitudes (or squared magnitudes) of the differences of these feature values is computed, and then added together using some user specified or automatically computed weights.

And although both simple and fast, there are two limitations to this notion of a "transition cost". The first problem is that we're mixing units (such as positions and velocities). This is why we need some kind of weighting terms to balance the contribution from each of these things. Setting good weights by hand is quite difficult, so we often resort to some kind of statistical normalization instead - something which can be a little brittle to changes in our data.

The second problem is that this does not account for the relationships between these features. And having some particular combination of joint position and velocity can often result in an unexpectedly good transition, which is not accounted for when we just add up those differences separately (as we will see later).

While I was thinking about this problem, I wondered if we could do something better by assuming our transitions are performed using inertialization. The idea I came up with was this: let's define our transition cost as the total displacement caused by an inertialized transition. i.e. The area of the graph between the source and destination animations:

inertialization example annotated

When we're using a critically damped spring to produce this offset between source and destination, something like this area can be computed directly using the integral of the spring equation.

For example, here's the integral of a critically damped spring, with target position zero and target velocity zero (as used to decay the offset of the inertializer), where the initial position is given as \( x \), the initial velocity as \( v \) and half-damping \( y \) is as defined here:

\begin{align*} \int (e^{-y \cdot t} \cdot (x + (v + x \cdot y) \cdot t))\, dt = -e^{-y \cdot t} \cdot \tfrac{1}{y^2} \cdot (t \cdot v \cdot y + x \cdot y \cdot (t \cdot y+2) + v)) \end{align*}

Then, by taking the definite integral of this equation from zero to infinity, we can get an equation which tells us the total area under the graph:

\begin{align*} \int_0^\infty (e^{-y \cdot t} \cdot (x + (v + x \cdot y) \cdot t))\, dt &= \\ -e^{-y \cdot \infty} \cdot \tfrac{1}{y^2} \cdot (\infty \cdot v \cdot y + x \cdot y \cdot (\infty \cdot y+2) + v)) &- \\ -e^{-y \cdot 0} \cdot \tfrac{1}{y^2} \cdot (0 \cdot v \cdot y + x \cdot y \cdot (0 \cdot y+2) + v)) &= \\ 0 + \tfrac{1}{y^2} \cdot (x \cdot y \cdot 2 + v)) &= \frac{2 \cdot x \cdot y + v}{y^2} \end{align*}

However this isn't exactly what we're looking for - and it wont exactly give us the area of the graph between source and destination animations in all cases. The first reason for this is that this number can be negative: springs which start with a negative displacement will produce a graph primarily below zero, resulting in a negative integral:

inertialization example offset 0

The second, worse problem, is that some graphs can even cross zero - producing an area both above and below zero, leading to an integral which underestimates the total area since it subtracts one from the other:

inertialization example offset 1

In these cases we need to split the integral into two parts - first computing the absolute value of the area before the graph crosses zero, and then the part after.

To find the time at which the graph crosses zero we can set the spring equation to equal zero and solve for time:

\begin{align*} e^{-y \cdot t} \cdot (x + (v + x \cdot y) \cdot t) = 0 \\ x + (v + x \cdot y) \cdot t = 0 \\ (v + x \cdot y) \cdot t = -x \\ t = \frac{-x}{v + x \cdot y} \end{align*}

The trick here is to remove the \( e^{-y \cdot t} \) term as for any value of \( t \) it simply scales the whole of the left hand side of the equation and so does not affect the point at which 0 is crossed.

Then, we need three different definite integrals, first from \( 0 \) to \( t \), from \( t \) to \( \infty \), and finally the one we already computed from \( 0 \) to \( \infty \):

\begin{align*} I(x, v, y, t) &= (e^{-y \cdot t} \cdot (x + (v + x \cdot y) \cdot t)) \\ \int_0^t I(x, v, y, t)\, dt &= \frac{2 \cdot x \cdot y + v}{y^2} - e^{-y \cdot t} \cdot \tfrac{1}{y^2} \cdot (t \cdot v \cdot y + x \cdot y \cdot (t \cdot y+2) + v)) \\ \int_t^\infty I(x, v, y, t)\, dt &= e^{-y \cdot t} \cdot \tfrac{1}{y^2} \cdot (t \cdot v \cdot y + x \cdot y \cdot (t \cdot y+2) + v)) \\ \int_0^\infty I(x, v, y, t)\, dt &= \frac{2 \cdot x \cdot y + v}{y^2} \\ \end{align*}

With these three integrals we can get an accurate estimate of the area of our transition. Given an initial positional difference x, velocity difference v and a halflife, we can compute the total displacement as follows:

def decay_spring_damper_intersection(x, v, halflife, eps=1e-8):
    y = halflife_to_damping(halflife) / 2.0
    return -x / (v + x * y)

def decay_spring_damper_displacement(x, v, halflife):
    y = halflife_to_damping(halflife) / 2.0
    t = decay_spring_damper_intersection(x, v, halflife)
    int_0 = (x*y*2 + v) / (y*y)
    int_t = (np.exp(-y*t) / (y*y)) * (t*v*y + x*y*(t*y+2) + v)
    return np.where(t > 0.0, 
        abs(int_0 - int_t) + abs(int_t),

But given we want to use this as a cost function, what does it actually look like? Well here's a 2D plot, with difference in position on the x-axis and difference in velocity on y-axis (the halflife here is set to 0.15):

inertialization cost function 2d

We can compare this to the cost function we would get if we just add the summed absolute differences in position and velocity (here the velocity difference is scaled by 0.5):

inertialization basic cost function 2d

The main difference we can see is that our new transition cost resembles more of a skewed valley than an inverted pyramid - the cost of a positive offset can actually be low - as long as it's combined with a negative velocity offset of the right magnitude. This makes sense - a large velocity offset can initialize the spring in a way such that it corrects itself back to zero more quickly than an initial velocity offset of zero:

inertialization example offset 2

But the problem with this cost function in practice is that since it's not a normal euclidean distance, it's not totally clear how we might integrate it with the rest of our system and exploit many of the common acceleration structures used in Motion Matching.

It would be nice if instead there were a "feature" we could compute which emulated the behavior of this cost function.

Well, as long as we're willing to accept the approximation \( \left| \tfrac{2 \cdot x \cdot y + v}{y^2} \right| \) which somewhat under-estimates the cost for oscillations around zero, there is!

The first step is to expand out the position difference \( x \) and velocity difference \( v \) from the source \( s \) and destination animations \( d \) into \( s_{pos} \), \( d_{pos} \), \( s_{vel} \), and \( d_{vel} \), giving:

\begin{align*} \left| \frac{2 \cdot (s_{pos} - d_{pos}) \cdot y + (s_{vel} - d_{vel})}{y^2} \right| \end{align*}

Then, we just need to re-arrange it until we've got all the \( s \) variables on one side of a difference, and all the \( d \) on the other:

\begin{align*} &= \left| \frac{2 \cdot y \cdot s_{pos} - 2 \cdot y \cdot d_{pos} + s_{vel} - d_{vel}}{y^2} \right| \\ &= \left| \frac{(2 \cdot y \cdot s_{pos} + s_{vel}) - (2 \cdot y \cdot d_{pos} + d_{vel})}{y^2} \right| \\ &= \left| \frac{(2 \cdot y \cdot s_{pos} + s_{vel})}{y^2} - \frac{(2 \cdot y \cdot d_{pos} + d_{vel})}{y^2} \right| \\ &= \left| \left( \frac{2 \cdot s_{pos}}{y} + \frac{s_{vel}}{y^2} \right) - \left( \frac{2 \cdot d_{pos}}{y} + \frac{d_{vel}}{y^2} \right) \right| \end{align*}

Which tells us that to emulate this approximate cost function the feature we need to put into our database is exactly the following:

\begin{align*} \frac{2 \cdot pos}{y} + \frac{vel}{y^2} \end{align*}

where \( pos \) is the bone position, \( vel \) is the bone velocity, and \( y \) is the half-damping which we compute from the half-life.

And although I'm certain there are some pathological cases which assign a small cost to transitions with large oscillations, in practice, if you use a relatively small half-life it seems to work. Here I've visualized the error of this function in comparison to our more exact cost function given previously.

inertialization approx cost function 2d

You can see that at least the error is limited to a fairly small part of the space where we have a large velocity offset and small, exact positional offset.

To test this feature, I added it to my previous Motion Matching demo, using the original setup and acceleration structure, but replacing the existing bone position and bone velocity features with this feature:

Not bad! And although the difference is not hugely significant, to me it looks like an overall improvement, and I was surprised at how well this idea worked out-of-the-box with no tweaking of velocity and positional weights required. So while perhaps not a silver bullet, to me this still seems like a fairly good way to potentially reduce the database size and remove one more weight for users to tweak!

As a disclaimer I should say that while I thought this was an interesting idea worth sharing, it's not something I've tested myself extensively, so in practice your mileage may vary.

Also worth noting is that although I've provided the derivation here for a spring-based inertializer, I think a similar derivation should be perfectly possible for inertializers that blend out the offset using a polynomial (and in fact it may well be easier).

Finally, it's worth mentioning that when we compute the total displacement here we are making an assumption that the inertialized movement of a joint's position in character space is going to be similar to what we would get if we individually inertialized all of the local joint rotations down the chain - something that may well not often be true.

Appendix: Cubic Inertializer

For fun I thought I would try and do the same derivation for a simple cubic inertializer function (which I used in my article on looping animations):

def decay_cubic(
    t = np.clip(dt / (blendtime + eps), 0, 1)

    d = x
    c = v * blendtime
    b = -3*d - 2*c
    a = 2*d + c
    return a*t*t*t + b*t*t + c*t + d

The idea behind this function is essentially just to find the coefficients of a cubic polynomial with some initial position and velocity on the y-axis, which crosses the x-axis at \( 1 \), with a velocity of \( 0 \). Then, if we want to change how long it takes to decay we can simply scale the x-axis (and adjust the initial velocity accordingly).

cubic intertialize function

The way we use this inertializer is a bit different. Rather than decaying the offset over time, we keep track of the initial offset, the time since transition, and compute the decayed offset on the fly based on how long it has been since the transition (we can actually do things this way with our spring-based inertializer too if we want).

void decayed_offset_cubic(
    vec3& out_x,
    vec3& out_v,
    const vec3 init_x,
    const vec3 init_v,
    const float blendtime, 
    const float dt,
    const float eps=1e-8f)
    float t = clampf(dt / (blendtime + eps), 0.0f, 1.0f);

    vec3 d = init_x;
    vec3 c = init_v * blendtime;
    vec3 b = -3.0f*d - 2.0f*c;
    vec3 a = 2.0f*d + c;
    out_x = a*t*t*t + b*t*t + c*t + d;
    out_v = (3.0f*a*t*t + 2.0f*b*t + c) / (blendtime + eps);

void inertialize_cubic_transition(
    vec3& off_x, 
    vec3& off_v, 
    float& off_t,
    const vec3 src_x,
    const vec3 src_v,
    const vec3 dst_x,
    const vec3 dst_v,
    const float blendtime)
    vec3 dec_x, dec_v;
    decayed_offset_cubic(dec_x, dec_v, off_x, off_v, blendtime, off_t);
    off_x = (dec_x + src_x) - dst_x;
    off_v = (dec_v + src_v) - dst_v;
    off_t = 0.0f;

void inertialize_cubic_update(
    vec3& out_x, 
    vec3& out_v,
    float& off_t,
    const vec3 off_x,
    const vec3 off_v,
    const vec3 in_x, 
    const vec3 in_v,
    const float blendtime,
    const float dt)
    off_t += dt;
    vec3 dec_x, dec_v;
    decayed_offset_cubic(dec_x, dec_v, off_x, off_v, blendtime, off_t);
    out_x = in_x + dec_x;
    out_v = in_v + dec_v;

Just like our spring, this function can cross the x-axis before \( 1 \), which means if we want to compute the absolute displacement, we need to compute the integral in two parts:

cubic intertialize intersection

And, just like before this means we need to compute the intersection time. In this case, since our cubic function is tangent to the x-axis at 1, we know that two of the roots will be exactly 1, which makes our lives a lot easier when it comes to finding the final one:

\begin{align*} a\ x^3 + b\ x^2 + c\ x + d = (x - 1)(x - 1)(x\ -\ ?) \end{align*}

To compute the final root we can use synthetic division:

\begin{align*} (a\ x^3 + b\ x^2 + c\ x + d)\ /\ (x - 1) &= a\ x^2 + (a + b)\ x + (a + b + c) \\ (a\ x^2 + (a + b)\ x + (a + b + c))\ /\ (x - 1) &= a\ x + 2a + b \end{align*}

Which means the value of the final root is given by the following:

\begin{align*} a\ x + 2a + b &= 0 \\ x &= \frac{-2\ a - b}{a} \end{align*}

We then just need to scale this by the blendtime if we want to get the actual intersection time:

def decay_cubic_intersection(x, v, blendtime, eps=1e-8):
    d = x
    c = v * blendtime
    b = -3*d - 2*c
    a = 2*d + c

    t = (-2*a - b) / (a + eps)
    return t * blendtime

So, depending on if there is an intersection between \( 0 \) and \( 1 \), the total displacement is either the blendtime multiplied by the integral of our polynomial between \( 0 \) and \( 1 \), or the blendtime multiplied by the sum of two integrals: one between \( 0 \) and \( t \), and one between \( t \) and \( 1 \) where \( t \) is the intersection time:

\begin{align*} \int_0^t a\ x^3 + b\ x^2 + c\ x + d\, dx &= \frac{a\ t^4}{4} + \frac{b\ t^3}{3} + \frac{c\ t^2}{2} + d\ t \\ \int_t^1 a\ x^3 + b\ x^2 + c\ x + d\, dx &= \left( \frac{a}{4} + \frac{b}{3} + \frac{c}{2} + d \right) - \left( \frac{a\ t^4}{4} + \frac{b\ t^3}{3} + \frac{c\ t^2}{2} + d\ t \right) \\ \int_0^1 a\ x^3 + b\ x^2 + c\ x + d\, dx &= \frac{a}{4} + \frac{b}{3} + \frac{c}{2} + d \\ \end{align*}

Which in code looks something like this:

def decay_cubic_displacement(x, v, blendtime, eps=1e-8):
    d = x
    c = v * blendtime
    b = -3*d - 2*c
    a = 2*d + c

    t = (-2*a - b) / (a + eps)
    int_0 = a/4 + b/3 + c/2 + d
    int_t = (a*t*t*t*t)/4 + (b*t*t*t)/3 + (c*t*t)/2 + (d*t)
    return blendtime * np.where((t >= 0.0) & (t <= 1.0),
        abs(int_0 - int_t) + abs(int_t),

Again, like before, if we can accept the approximation which does not split the integral into two parts, and underestimates the displacement of oscillations, we can derive a motion matching feature by getting all the source variables and destinations variables on either side of a subtraction, where here \( z \) is the blendtime:

\begin{align*} &= z\ \left| \frac{a}{4} + \frac{b}{3} + \frac{c}{2} + d \right| \\ &= z\ \left| \frac{z}{2}\ s_{vel} - \frac{z}{2}\ d_{vel} - \frac{2 z}{3}\ s_{vel} + \frac{2 z}{3}\ d_{vel} + \frac{1}{2}\ s_{pos} - \frac{1}{2}\ d_{pos} + \frac{z}{4}\ s_{vel} - \frac{z}{4}\ d_{vel} \right| \\ &= z\ \left| \left( \frac{z}{2}\ s_{vel} - \frac{2 z}{3}\ s_{vel} + \frac{1}{2}\ s_{pos} + \frac{z}{4}\ s_{vel} \right) - \left( \frac{z}{2}\ d_{vel} - \frac{2 z}{3}\ d_{vel} + \frac{1}{2}\ d_{pos} + \frac{z}{4}\ d_{vel} \right) \right| \\ &= z\ \left| \left( \frac{z}{12}\ s_{vel} + \frac{1}{2}\ s_{pos} \right) - \left( \frac{z}{12}\ d_{vel} + \frac{1}{2}\ d_{pos} \right) \right| \\ \end{align*}

Giving the following feature:

\begin{align*} \frac{z \cdot pos}{2} + \frac{z^2\cdot vel}{12} \end{align*}

where \( pos \) is the bone position, \( vel \) is the bone velocity, and \( z \) is the blendtime. I find it interesting how close this is to the spring-based feature given that the half-damping \( y \) is kind of similar to an inverse blend time.

Finally, here you can see an adaption of my motion matching demo, with this feature used for matching, and the cubic inertializer used for blending:


github twitter rss