Cubic Interpolation of Quaternions

2023-05-03 01:27:05

Created on April 25, 2023, 9:01 a.m.

If you happen to’ve ever googled “cubic interpolation of quaternions” or appeared up the “SQUAD” algorithm you would be forgiven for pondering that understanding methods to do clean, cubic interpolation of quaternions requires some actually superior arithmetic.

But it surely actually would not have to be that sophisticated. And if we perceive properly the algorithm for regular cubic interpolation (and extra particularly, Catmull-Rom cubic interpolation), the quaternion model just about falls into our laps.

Each time the topic of cubic interpolation and splines comes up it might be incorrect of me to not suggest absolutely the master-class that’s Freya Holmér’s spline videos – watching these ought to provide you with an ideal instinct not simply in regards to the course of we’re going to undergo on this article, however about splines normally.

Hermite Curve Interpolation

Neglect about quaternions for now. As an alternative, we’ll begin with the next extra easy downside: if now we have two values at 0 and 1, and gradients for these values at 0 and 1, how can we produce a clean interpolation which passes via these values with the given gradients?

The answer is to suit a cubic polynomial perform ( f ) to the 2 factors and their gradients. On this case now we have 4 constraints…

start{align*}
f(0) &= p_0
f(1) &= p_1
f'(0) &= v_0
f'(1) &= v_1
finish{align*}

…and 4 unknowns ( a ), ( b ), ( c ), and ( d ).

start{align*}
f(x) &= a x^3 + b x^2 + c x + d
f'(x) &= 3 a x^2 + 2 b x + c
finish{align*}

So we simply must plug our totally different equations into one another and resolve for ( a ), ( b ), ( c ), and ( d ).

First we resolve for ( d ) utilizing our first constraint ( f(0) = p_0 )…

start{align*}
p_0 &= a 0 + b 0 + c 0 + d
d &= p_0
finish{align*}

Then we resolve for ( c ) utilizing our third constraint ( f'(0) = v_0 )…

start{align*}
v_0 &= 3 a 0 + 2 b 0 + c
c &= v_0
finish{align*}

Then we will resolve for ( b ) utilizing our different two constraints and our options for ( d ) and ( c )…

start{align*}
p_1 &= a 1 + b 1 + c 1 + d
a &= p_1 – b – v_0 – p_0
&
v_1 &= 3 a 1 + 2 b 1 + c
v_1 &= 3 (p_1 – b – v_0 – p_0) + 2 b + v_0
v_1 &= 3 p_1 – 3 b – 3 v_0 – 3 p_0 + 2 b + v_0
b &= 3 p_1 – 2 v_0 – 3 p_0 – v_1
finish{align*}

And eventually we will resolve for ( a )…

start{align*}
a &= p_1 – b – v_0 – p_0
a &= p_1 – (3 p_1 – 2 v_0 – 3 p_0 – v_1) – v_0 – p_0
a &= v_0 + 2 p_0 + v_1 – 2 p_1
finish{align*}

We are able to plot this cubic perform to confirm that it passes via the given factors, and with the right gradients:

Good! If we broaden the cubic perform with our values for ( a ), ( b ), ( c ), and ( d ) substituted, we will rearrange issues such that our result’s expressed when it comes to ( p_0 ), ( p_1 ), ( v_0 ) and ( v_1 ) instantly (quite than the polynomial coefficients):

start{align*}
f(x) &= (v_0 + 2 p_0 + v_1 – 2 p_1) x^3 + (3 p_1 – 2 v_0 – 3 p_0 – v_1) x^2 + v_0 x + p_0
f(x) &= v_0 x^3 + 2 p_0 x^3 + v_1 x^3 – 2 p_1 x^3 + 3 p_1 x^2 – 2 v_0 x^2 – 3 p_0 x^2 – v_1 x^2 + v_0 x + p_0
f(x) &= p_0 (2 x^3 – 3 x^2 + 1) + p_1 (3 x^2 – 2 x^3) + v_0 (x^3 – 2 x^2 + x) + v_1 (x^3 – x^2)
finish{align*}

That is helpful as a result of often ( x ) is a scalar, whereas ( p_0 ), ( p_1 ), ( v_0 ) and ( v_1 ) can usually be 3d vectors. For instance, we will implement a model that works on 3d vectors in C++ as follows, the place x is our interpolating worth between 0 and 1:

vec3 hermite_basic(float x, vec3 p0, vec3 p1, vec3 v0, vec3 v1)
{  
    float w0 = 2*x*x*x - 3*x*x + 1;
    float w1 = 3*x*x - 2*x*x*x;
    float w2 = x*x*x - 2*x*x + x;
    float w3 = x*x*x - x*x;
    
    return w0*p0 + w1*p1 + w2*v0 + w3*v1;
}

(Apart: Discover how w0 and w1 are each variations of smoothstep – is not that cool!)

Since now we have the spinoff of our cubic polynomial perform, we will additionally compute the gradient/velocity of the interpolated end result and get our perform to output that too:

void hermite(
    vec3& pos,
    vec3& vel, 
    float x, 
    vec3 p0,
    vec3 p1, 
    vec3 v0,
    vec3 v1)
{
    float w0 = 2*x*x*x - 3*x*x + 1;
    float w1 = 3*x*x - 2*x*x*x;
    float w2 = x*x*x - 2*x*x + x;
    float w3 = x*x*x - x*x;
    
    float q0 = 6*x*x - 6*x;
    float q1 = 6*x - 6*x*x;
    float q2 = 3*x*x - 4*x + 1;
    float q3 = 3*x*x - 2*x;
    
    pos = w0*p0 + w1*p1 + w2*v0 + w3*v1;
    vel = q0*p0 + q1*p1 + q2*v0 + q3*v1;
}

Which appears one thing like this:

One final adjustment to this formulation is to successfully translate our factors such that ( p_0 = 0 ). Then, as soon as now we have the interpolated worth, translate it again up by including the unique worth of ( p_0 ):

void hermite_alt(
    vec3& pos,
    vec3& vel, 
    float x, 
    vec3 p0,
    vec3 p1, 
    vec3 v0,
    vec3 v1)
{
    vec3 p1_sub_p0 = p1 - p0;

    float w1 = 3*x*x - 2*x*x*x;
    float w2 = x*x*x - 2*x*x + x;
    float w3 = x*x*x - x*x;
    
    float q1 = 6*x - 6*x*x;
    float q2 = 3*x*x - 4*x + 1;
    float q3 = 3*x*x - 2*x;
    
    pos = w1*p1_sub_p0 + w2*v0 + w3*v1 + p0;
    vel = q1*p1_sub_p0 + q2*v0 + q3*v1;
}

This may seem to be a pointless further step, however this formulation with the primary level at zero really simplifies issues in a approach which will likely be actually helpful once we begin coping with quaternions.

What we have derived is known as Hermite Cubic Interpolation, and is actually the arithmetic used for a single phase of a Hermite Spline, which will likely be a key a part of doing cubic interpolation of quaternions.

Catmull-Rom Cubic Interpolation

If now we have a sequence of factors and their gradients/velocities, we will use the earlier equation to make a clean interpolation that passes via them. However what if we simply have factors? With no velocities or gradients related to them? Can we nonetheless use Hermite Curve Interpolation?

Sure! The trick is to take 4 factors as a substitute of two, and to make use of the central difference to “estimate” a gradient/velocity to affiliate with every of the center two factors. Then, as soon as now we have these velocities we will use Hermite Curve Interpolation to interpolate the intermediate values.

In code it appears one thing like this:

void catmull_rom(
    vec3& pos,
    vec3& vel,
    float x,
    vec3 p0,
    vec3 p1, 
    vec3 p2,
    vec3 p3)
{
    vec3 v1 = ((p1 - p0) + (p2 - p1)) / 2;
    vec3 v2 = ((p2 - p1) + (p3 - p2)) / 2;
    return hermite(pos, vel, x, p1, p2, v1, v2);
}

You’ll be able to see that we compute the central distinction for these two center factors utilizing the typical of the ahead and backward variations. Then, we move these estimated velocities to our Hermite Curve Interpolation perform to get the ultimate interpolated worth! The result’s one thing like this:

And that is how we use a Catmull-Rom spline to do cubic interpolation. Here’s what it appears like in 3d:

Quaternion Catmull-Rom Cubic Interpolation

Now that we perceive precisely what’s going on (and what every totally different worth and computation represents) we’re able to adapt our equations to work with quaternions as a substitute of positions.

First we’ll sort out the Hermite Cubic Interpolation perform, and there are a couple of modifications we have to make to get this formulation to work for quaternions:

As an alternative of linear velocities, we have to use angular velocities.

Once we compute the distinction between p0 and p1 to place p0 at zero, we’re as a substitute going to take the quaternion distinction.

As soon as now we have this quaterion distinction, we’re going to convert it into the scaled-angle-axis area in order that it is sensible to combine it with angular velocities.

As soon as now we have added every part collectively, we have to convert the end result again to a quaternion. And as a substitute of including again the unique worth of p0 we’ll use quaternion multiplication.

In code it appears like this:

void quat_hermite(
    quat& rot,
    vec3& vel, 
    float x, 
    quat r0,
    quat r1, 
    vec3 v0,
    vec3 v1)
{
    float w1 = 3*x*x - 2*x*x*x;
    float w2 = x*x*x - 2*x*x + x;
    float w3 = x*x*x - x*x;
    
    float q1 = 6*x - 6*x*x;
    float q2 = 3*x*x - 4*x + 1;
    float q3 = 3*x*x - 2*x;
    
    vec3 r1_sub_r0 = quat_to_scaled_angle_axis(quat_abs(quat_mul_inv(r1, r0)));   
    
    rot = quat_mul(quat_from_scaled_angle_axis(w1*r1_sub_r0 + w2*v0 + w3*v1), r0);
    vel = q1*r1_sub_r0 + q2*v0 + q3*v1;
}

Subsequent, we will sort out our Catmull-Rom Cubic Interpolation perform. This time there are fewer changes. All we actually must do is convert our velocity computations through finite distinction, into angular velocity computations through finite distinction (as described in my previous article on angular velocity)…

void quat_catmull_rom(
    quat& rot,
    vec3& vel,
    float x,
    quat r0,
    quat r1, 
    quat r2,
    quat r3)
{
    vec3 r1_sub_r0 = quat_to_scaled_angle_axis(quat_abs(quat_mul_inv(r1, r0)));
    vec3 r2_sub_r1 = quat_to_scaled_angle_axis(quat_abs(quat_mul_inv(r2, r1)));
    vec3 r3_sub_r2 = quat_to_scaled_angle_axis(quat_abs(quat_mul_inv(r3, r2)));
  
    vec3 v1 = (r1_sub_r0 + r2_sub_r1) / 2;
    vec3 v2 = (r2_sub_r1 + r3_sub_r2) / 2;
    return quat_hermite(rot, vel, x, r1, r2, v1, v2);
}

And that is it! Whereas a few of these steps won’t be 100% apparent with out enthusiastic about it somewhat, so long as you perceive the connection between angular velocities and quaternions, I do not assume there may be something surprising happening.

Here’s what it appears like in our little 3d atmosphere:

Uncooked Quaternion Cubic Interpolation

If you’re interpolating a sequence of quaternions that are sequentially fairly related to one another, and have been “unrolled” in order that there are not any sudden
discontinuities there may be a fair simpler option to do cubic interpolation of quaternions, which is to simply do the interpolation instantly within the 4d quaternion area, treating the quaternions as 4d vectors, and re-normalize the end result.

This may be quicker and can produce virtually similar outcomes when your quaternions are related. I see this as like the difference between using slerp and nlerp for linear interpolation.

Conclusion

In the case of animation knowledge, I believe most individuals consider that doing cubic interpolation of rotations is a waste of CPU cycles. Sure – sampling 4 poses is dearer than two and the arithmetic concerned are dearer too – however cubic interpolation actually can make a visible distinction – specifically when taking part in animation in slow-motion.

As well as, in comparison with linear interpolation, cubic interpolation offers velocities that change easily in-between frames, which might forestall aliasing results when additional processing the info. For instance, sampling velocities for a motion-matching database through linear interpolation at a fee increased than the unique animation knowledge will produce consecutive entries with the identical velocities – which might seem like a bug in the case of examine the info and might have an effect on the results of downstream algorithms akin to PCA.

Right here is one other concrete instance: whereas FIFA 21 has undoubtedly a few of the finest, most refined, and most life like animation on the earth of video video games, the linear interpolation between frames, and the speed discontinuity this introduces throughout slow-mo is one thing that can’t be un-seen.