Icon Quaternion Weighted Average

 

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 RAGE

 

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 LOL I DREW THIS DRAGON

 

Icon One Cat Just Leads To Another

Quaternion Weighted Average

Created on Dec. 31, 2023, 10:52 a.m.

How do we find the weighted average of a set of quaternions?

This is not immediately obvious, in particular given that a quaternion and its negation represent the same rotation - so just doing a normal weighted sum of the raw 4d quaternion values is not going to work.

Well, here is a simple solution you might have seen used in various places:

quat quat_average_basic(
    const slice1d<quat> rotations, 
    const slice1d<float> weights)
{
    quat accum = quat(0.0f, 0.0f, 0.0f, 0.0f);
    
    // Loop over rotations
    for (int i = 0; i < rotations.size; i++)
    {
        // If more than 180 degrees away from current running average...
        if (quat_dot(accum, rotations(i)) < 0.0f)
        {
            // Add opposite-hemisphere version of the quaternion
            accum = accum - weights(i) * rotations(i);
        }
        else
        {
            // Add the quaternion
            accum = accum + weights(i) * rotations(i);
        }
    }
    
    // Normalize the final result
    return quat_normalize(quat_abs(accum));
}

Essentially we just keep a running average and flip any quaternions that are on the opposite hemisphere to our running average. We then normalize the final result, and can take the version of this average rotation closest to the identity rotation using the quat_abs function if we like.

This is simple and works pretty well in most cases - in particular if the rotations in the set are already similar.

There are, however, a few problems we might encounter with this method: first - because we are using a running average the result of this function can actually change depending on the order of the provided quaternions.

For example, here is the rotation we get (shown on the right) if we apply this averaging function to sets of randomly generated quaternions (shown on the left) with a toggle to reverse the order in which we provide the quaternions to the function:

As you can see, sometimes the result is the same, but sometimes it can be really quite different. Not a desirable property of an averaging function!

Another issue with this method is that the result it produces has many unnecessary discontinuities, and can often suddenly jump to a new result as we change the rotations in the set. We can see this if we slowly adjust one of the rotations in the set:

For large sets of quaternions that are similar to each other this might not be an issue. But for smaller sets of quaternions or quaternions with large angular differences this can (like the nlerp function) make sudden jumps when the hemisphere of one of the quaternions changes. In fact, if we run this function on just two quaternions, the result is identical to nlerp.

So is there a better way to compute the average rotation? Yes - although first we need to think a bit more about what we actually want an average function to give us. One idea is this: we want the average rotation of a set of rotations to be a rotation which is as close as possible to all the rotations in that set. More specifically we can say that we want the average rotation of a set to be a rotation which has the smallest sum of squared distances to all the rotations in that set.

The derivation for how to compute this is a little complex (described in this paper), but the final implementation can be relatively simple. First we compute a weighted sum of the outer product of all our quaternions as a 4x4 matrix:

mat4 accum = mat4_zero();

// Loop over rotations
for (int i = 0; i < rotations.size; i++)
{
    quat q = rotations(i);
     
    // Compute the outer-product of the quaternion 
    // multiplied by the weight and add to the accumulator 
    accum = accum + weights(i) * mat4(
        q.w*q.w, q.w*q.x, q.w*q.y, q.w*q.z,
        q.x*q.w, q.x*q.x, q.x*q.y, q.x*q.z,
        q.y*q.w, q.y*q.x, q.y*q.y, q.y*q.z,
        q.z*q.w, q.z*q.x, q.z*q.y, q.z*q.z);
}

Once we have this, the average quaternion as given by the above definition, will be exactly equal to the largest eigenvector of this 4x4 matrix. This can be computed using the power iteration method (which you may remember from my previous article on Joint Limits):

vec4 mat4_svd_dominant_eigen(
    const mat4 A, 
    const vec4 v0,
    const int iterations, 
    const float eps)
{
    // Initial Guess at Eigen Vector & Value
    vec4 v = v0;
    float ev = (mat4_mul_vec4(A, v) / v).x;
    
    for (int i = 0; i < iterations; i++)
    {
        // Power Iteration
        vec4 Av = mat4_mul_vec4(A, v);
        
        // Next Guess at Eigen Vector & Value
        vec4 v_new = normalize(Av);
        float ev_new = (mat4_mul_vec4(A, v_new) / v_new).x;
        
        // Break if converged
        if (fabs(ev - ev_new) < eps)
        {
            break;
        }
        
        // Update best guess
        v = v_new;
        ev = ev_new;
    }
    
    return v;
}

This gives the following more accurate quaternion averaging function:

quat quat_average_accurate(
    const slice1d<quat> rotations, 
    const slice1d<float> weights)
{
    mat4 accum = mat4_zero();
    
    // Loop over rotations
    for (int i = 0; i < rotations.size; i++)
    {
        quat q = rotations(i);
      
        // Compute the outer-product of the quaternion 
        // multiplied by the weight and add to the accumulator 
        accum = accum + weights(i) * mat4(
            q.w*q.w, q.w*q.x, q.w*q.y, q.w*q.z,
            q.x*q.w, q.x*q.x, q.x*q.y, q.x*q.z,
            q.y*q.w, q.y*q.x, q.y*q.y, q.y*q.z,
            q.z*q.w, q.z*q.x, q.z*q.y, q.z*q.z);
    }
    
    // Initial guess at eigen vector is identity quaternion
    vec4 guess = vec4(1, 0, 0, 0);
    
    // Compute first eigen vector
    vec4 u = mat4_svd_dominant_eigen(accum, guess, 64, 1e-5f);
    vec4 v = normalize(mat4_transpose_mul_vec4(accum, u));
    
    // Average quaternion is first eigen vector
    return quat_abs(quat(v.x, v.y, v.z, v.w));
}

The distance this algorithm minimizes is effectively the squared distance between rotation matrices. i.e. it minimizes the sum of the squared distances between the rotation matrix of this average quaternion, and the rotation matrices of every quaternion in the set. And since we are talking in terms of distances between matrices, specifically it minimizes the sum of the squared Frobenius distances.

For rotations, the Frobenius distance between the rotation matrices is actually equivalent to the sin of half the angle between the two rotations:

float quat_frobenius_distance(quat q, quat p)
{
    // Compute the difference between two quaternions
    quat d = quat_abs(quat_mul_inv(q, p));
    
    // Compute the half-angle of that difference
    float halfangle = acosf(clampf(d.w, -1.0f, 1.0f));
    
    // Compute the sin of the half angle
    return sinf(halfangle);
}

This squared sin of the half-angle may seem like an odd kind of distance to be minimizing. It might seem more intuitive to try to minimize the sum of the squared angular differences directly, however the square of the sin of the half angle has some nice properties - it is similar to the square of the half-angle itself for small angles, but it plateaus smoothly around an angular difference of pi radians (i.e. when the two rotations are at their most different), which gives us smooth continuous behavior for the averaging function:

quat average distance

Although more expensive to compute (and the constant time overhead is quite large if you have just a few rotations in the set), this accurate algorithm for averaging quaternions solves our issues with the previous method. It's independent of the order we provide the rotations:

And it does not jump when we slowly adjust the rotations:

To test that this algorithm does what it states, we can compute the sum of the squared sin-half-angles between our average rotation and the rotations in the set. This will help us verify that this algorithm really does produce a smaller result than our basic method:

So yes, this method really does produce an average rotation with a smaller sum of squared sin-half-angles than our first basic method.

If you want to play around with the demo shown in this article the source code can be found here.

And that's about it. If you want to read more, there are lots of other resources online covering the same topic from which most of this article is taken! And as usual, thanks for reading.


Appendix: Weighted Average of Relative Rotations

As I discussed in my article on Joint Limits, it's really important when dealing with quaternions to keep in mind if they represent relative rotations or absolute rotations - i.e. are these rotations of something relative to something else meaningful, or are they just orientations in the world space that are not relative to anything in particular.

And as also discussed in that article - if we have joint rotations representing a character's pose, what these rotations are relative to is the pose we get when we set all rotations to the identity: the identity pose - something that is usually not at all a meaningful pose (which makes these rotations more like absolute rotations):

Just as how applying joint limits with respect to the identity pose does not work - just taking the raw average of quaternions with respect to the identity pose does not work either - and is what will produce the discontinuities and bad results we saw before.

void bone_rotations_weighted_average_raw(
    slice1d<quat> blended_rotations, 
    const slice2d<quat> sample_rotations, 
    const slice1d<float> sample_weights)
{
    assert(sample_rotations.rows == sample_weights.size);
    assert(sample_rotations.cols == blended_rotations.size);
    
    blended_rotations.zero();
    
    for (int i = 0; i < sample_rotations.rows; i++)
    {
        for (int j = 0; j < sample_rotations.cols; j++)
        {
            blended_rotations(j) = blended_rotations(j) + 
                sample_weights(i) * quat_abs(sample_rotations(i, j));
        }
    }
    
    for (int j = 0; j < sample_rotations.cols; j++)
    {
        blended_rotations(j) = quat_normalize(blended_rotations(j));
    }
}

For example, this is the result we get if we take the raw weighted average of the quaternion values of these three different poses, as compared to the "accurate" method I presented above:

Similarly, we could do the weighted average in the log space instead...

void bone_rotations_weighted_average_log(
    slice1d<quat> blended_rotations, 
    slice1d<vec3> accum_rotations, 
    const slice2d<quat> sample_rotations, 
    const slice1d<float> sample_weights)
{
    assert(sample_rotations.rows == sample_weights.size);
    assert(sample_rotations.cols == blended_rotations.size);
    assert(sample_rotations.cols == accum_rotations.size);
    
    accum_rotations.zero();
    
    for (int i = 0; i < sample_rotations.rows; i++)
    {
        for (int j = 0; j < sample_rotations.cols; j++)
        {
            accum_rotations(j) += sample_weights(i) * 
                quat_log(quat_abs(sample_rotations(i, j)));
        }
    }
    
    for (int j = 0; j < sample_rotations.cols; j++)
    {
        blended_rotations(j) = quat_exp(accum_rotations(j));
    }
}

but this is not any better...

The problem is that these rotations are being averaged "around" this identity pose - a pose which is not meaningful to take an average around - and just moving to the log space does not change this.

If we instead transform all our rotations to be relative to some meaningful character pose such as the skinned rest pose...

...then taking the weighted average makes a lot more sense and gives a lot better results.

void bone_rotations_weighted_average_raw_ref(
    slice1d<quat> blended_rotations, 
    const slice1d<quat> reference_rotations,
    const slice2d<quat> sample_rotations, 
    const slice1d<float> sample_weights)
{
    assert(sample_rotations.rows == sample_weights.size);
    assert(sample_rotations.cols == blended_rotations.size);
    assert(sample_rotations.cols == reference_rotations.size);
    
    blended_rotations.zero();
    
    for (int i = 0; i < sample_rotations.rows; i++)
    {
        for (int j = 0; j < sample_rotations.cols; j++)
        {
            blended_rotations(j) = blended_rotations(j) + sample_weights(i) * 
                quat_abs(quat_inv_mul(
                    reference_rotations(j), sample_rotations(i, j)));
        }
    }
    
    for (int j = 0; j < sample_rotations.cols; j++)
    {
        blended_rotations(j) = quat_abs(quat_mul(
            reference_rotations(j), quat_normalize(blended_rotations(j))));
    }
}

Here is what happens visually if we make these rotations relative to the "reference pose" before taking the average:

Now we can see that the result is really similar to our "accurate" method from before.

This method is faster to compute since we don't need to compute an eigenvector which is very nice. Also, in some cases it might even give more desirable results since it will effectively average a rotation of -170 degrees away from the reference pose in one direction, and 170 degrees away from the reference pose in the opposite direction back to the reference pose itself. Our "accurate" method on the other hand would average these two rotations to a full rotation of 180 degrees. Which behavior you prefer here obviously depends on the application.

Notice that the result of our "accurate" method is unchanged by using rotations relative to the reference pose or not - it always treats rotations as absolute rotations, exactly as hoped for!

Due to the way the quaternion space is shaped, the weighted average in the raw quaternion space will tend to reduce the impact of large rotations close to 180 degrees. If we want to take these into account more fairly it is better to do the averaging in the log space.

void bone_rotations_weighted_average_log_ref(
    slice1d<quat> blended_rotations, 
    slice1d<vec3> accum_rotations, 
    const slice1d<quat> reference_rotations,
    const slice2d<quat> sample_rotations, 
    const slice1d<float> sample_weights)
{
    assert(sample_rotations.rows == sample_weights.size);
    assert(sample_rotations.cols == blended_rotations.size);
    assert(sample_rotations.cols == accum_rotations.size);
    assert(sample_rotations.cols == reference_rotations.size);
    
    accum_rotations.zero();
    
    for (int i = 0; i < sample_rotations.rows; i++)
    {
        for (int j = 0; j < sample_rotations.cols; j++)
        {
            accum_rotations(j) += sample_weights(i) * 
                quat_log(quat_abs(quat_inv_mul(
                    reference_rotations(j), sample_rotations(i, j))));
        }
    }
    
    for (int j = 0; j < sample_rotations.cols; j++)
    {
        blended_rotations(j) = quat_abs(quat_mul(
            reference_rotations(j), quat_exp(accum_rotations(j))));
    }
}

However, the difference is fairly small in this example - and the computational cost is slightly larger:

If you're interested in learning more about making rotations relative to a reference pose, check out this old video by Casey Muratori.

To summarize:

If you have a "reference" you want to average around then transform your rotations relative to that reference and either average the raw quaternion values, or average in the log space (if you want slightly better quality at the cost of some performance loss).

If you have no "reference" or are averaging quaternions that are not relative to anything meaningful, use the "accurate" quaternion averaging function as presented above.

In fact, if you have no "reference" but know the quaternions in your set are going to be similar to some other set, one approach I've seen before is to use the "accurate" method ahead of time on this other set to extract a reference which future rotations can be more performantly averaged around later on. So that's an option too!

github twitter rss