Maths - Quaternion Interpolation (SLERP) - Forum discussion

By: martinbaker ( Martin Baker )
quaternion spherical interpolation SLERP  
2003-09-05 17:44

Following on from the interesting ideas from minorlogic on the 'quaternion from 2 vec' thread, I have added a new page:

http://www.euclideanspace.com/maths/algebra/realNormedAlgebra/quaternions/slerp/

I would appreciate any help in developing the page (it has some questions in it).

Thanks,

Martin

By: s_ludwig ( Sven Ludwig )
RE: quaternion spherical interpolation SLERP  
2003-09-06 08:59

Hi,

> cos(theta) = qa.x * qb.x + qa.y * qb.y + qa.z * qb.z + qa.w * qb.w
> Do you know where this formula comes from?

if I am not mistaken, if quaternions are represented
by a scalar w and a three-vector v = [x,y,z]^T as (w,v),
then this formula is the scalar w of the product of
qa and the conjugate of qb, or vice-versa the product
of the conjugate of qa with qb.

Cheers,
Sven

By: martinbaker ( Martin Baker )
RE: quaternion spherical interpolation SLERP  
2003-09-06 15:32

Sven,

I think I see what you mean, if q is the quaternion which rotates from qa to qb then:
qb = qa * q

multiplying both sides by conj(qa) gives:
q = conj(qa) * qb

The real part of a multiplication is:
real((qa.w + i qa.x + j qa.y + k qa.z)*(qb.w + i qb.x + j qb.y + k qb.z)) = qa.w * qb.w - qa.x*qb.x - qa.y*qb.y- qa.z*qb.z

So using the conjugate of qa gives:
real((qa.w - i qa.x - j qa.y - k qa.z)*(qb.w + i qb.x + j qb.y + k qb.z)) = qa.w*qb.w + qa.x*qb.x + qa.y*qb.y+ qa.z*qb.z

real(q) = cos(t/2)

Therefore
cos(t/2) = qa.w*qb.w + qa.x*qb.x + qa.y*qb.y+ qa.z*qb.z

Thanks,

Martin

By: minorlogic ( Michaele Norel )
RE: quaternion spherical interpolation SLERP  
2003-09-09 04:18

Hi Martin!
Can't find the Shoemake paper too , but remember that sow many provs of SLERP.

I try (some time) to find the more efficient way to do SLERP and avoid some bad math cases. Doing this found interest page:
http://www.hadron.org/~hatch/rightway.php

By: minorlogic ( Michaele Norel )
RE: quaternion spherical interpolation SLERP  
2003-09-09 04:32

But Don Hatch miss something:
The thing that he try to do more robust if the :

( q1*sin((1-t)*theta) + q2*sin(t*theta))/sin(theta) can be just scaled by sin(theta) and than result normalized.

( q1*sin((1-t)*theta) + q2*sin(t*theta))

May be you will find some less trigonometry ? To be stable near 0 angle i propose to multiply quaternion by left or right part :

( q1*sin((1-t)*theta)/sin(t*theta) + q2) and than make some trig manipulation to make it easier.

q1(sin(Theta)*ctg( Theta*t ) - cos(Theta)) + q2 - for example

Maybe somebody can help ?

By: martinbaker ( Martin Baker )
RE: quaternion spherical interpolation SLERP  
2003-09-09 10:00

Hi minorlogic,

Yes the rightway website is very good, I would like to go over all the algorithms on my website and check them for stability in that level of detail. I think that's a long term project though!

Do you think there is a way to use the fact that we are starting with cos(theta) rather than theta? We can go straight to sin(theta) by using Pythagoras s = sqrt(1 - c^2). But is there a way to get to sin(t*theta) or sin((1-t)*theta) without calculating theta first?

Martin

By: minorlogic ( Michaele Norel )
RE: quaternion spherical interpolation SLERP  
2003-09-10 03:41

Hi Martin !
Yes i hope the cos(theta) can be used in more efficient way.

I don't know how to find sin((1-t)*theta) or ctg( theta*t ), i think about , may be exist the other way. May be it fundamental problem ? We interpolate the angels but have only thair projections ( sin , cos ) and we can't solve it without projection->angle->projection conversion ?
// ----
Exist a very efficient approximation to this, one of them was published :
the article "Hacking Quaternions",
The Inner Product, by Jonathan Blow, Game Developer Magazine,
March 2002.
And source code can be downloaded.

By: minorlogic ( Michaele Norel )
RE: quaternion spherical interpolation SLERP  
2003-10-09 03:34

Hi Martin !

Thinking about the "The Right Way to Calculate Stuff" Article? and the way to calculate the angle between vectors.

angle = atan2( |v1Xv2|, v1.v2 )

found that for example the quaternion -> axis angle can be done same way

angle = atan2( |quat.vec|, quat.w )
(and quat can be nonunit )

And the same can be done with SLERP to calculate the angle between quats but stable.

float calculate_theta(quaternion& q1, quaternion& q2 )
{
quaternion tq = q1;
tq.conjugate();
tq = tq * q2;

return atan2( sqrt(tq.x*tq.x + tq.y*tq.y + tq.z*tq.z), tq.w );
}
I think this is the way to get stable calculations.

What do you think about ?

By: martinbaker ( Martin Baker )
RE: quaternion spherical interpolation SLERP  
2003-10-09 09:57

Minorlogic,

Thanks, is there a proof of: angle = atan2( |quat.vec|, quat.w ) ?
And is it more stable than: angle = 2 * acos(w) ?

I can almost see a proof, as follows, but it doesn't quite work? The divide by 2 disappears at the last stage, I think it needs a trig identity for tan(angle/2)?

tan(angle/2) = sin(angle/2)/cos(angle/2) // ie trig identity tan = sin/cos

but |quat.vec| = sin(angle/2)
and w = cos(angle/2)

which gives:
tan(angle/2) =|quat.vec| / cos(angle/2)
angle = atan2( |quat.vec|, quat.w )

Martin

By: minorlogic ( Michaele Norel )
RE: quaternion spherical interpolation SLERP  
2003-10-10 02:11

Thanks Martin.
Yes my mistake !
The rotation angle :

angle = 2.0 * atan2( |quat.vec|, quat.w ) - for axis angle conversion

But i thinked about SLERP, when writed this , and there is used HALF angle

half_angle = atan2( |quat.vec|, quat.w )

tan(a/2) = sin(a/2)/cos(a/2)

if quat is not unit the scale goes when dividing.

q( s * sin(a/2)*v, s * cos(a/2) )
(s*sin(a/2))/(s*cos(a/2)) = sin(a/2)/cos(a/2)

//-------------
"Thanks, is there a proof of: angle = atan2( |quat.vec|, quat.w ) ?
And is it more stable than: angle = 2 * acos(w) ?"

Not sure about proof, but easy can see that near

cos -> 1 or cos -> -1 the precision of acos goes to nothing. (the derivative goes to infiniti) , im not sure that proof of it is easy...

I use intuition and experemental data. May be you can proofe it ?

By: minorlogic ( Michaele Norel )
RE: quaternion spherical interpolation SLERP  
2003-10-22 03:28

Hi Martin.
I made some tests, and found that:
function when calculating angle between 2 vectors
version that use atan2 3 time faster than version that use acos.

And experimental error with "atan2" always less than with "acos"

Some peoples pointed me that atan2 is accelerated by x86 asm "fpatan", but acos not. And the acos calculated by some:

atan2( sqrt(1 - val*val), val );

May be it is intresting.

By: martinbaker ( Martin Baker )
RE: quaternion spherical interpolation SLERP  
2003-10-22 08:07

Hi minorlogic,

Yes that's interesting, my first concern was what happens at the singularity of the tan function at 90 degrees when tan(90) jumps from +infinity to -infinity, but thinking about it I guess that's not a problem for atan function (only tan function).

More of a problem might be 0 degrees where both tan(small number) and cos(small number) are flat so atan(near zero) and acos(near +-one) will produce the greatest fluctuations in angle and therefore be least accurate? Does this agree with your tests? I guess this is not really an issue here, even if both functions are least accurate at zero degrees atan will still be better?

Martin.

By: minorlogic ( Michaele Norel )
RE: quaternion spherical interpolation SLERP  
2003-10-22 08:47

Hi Martin.
The all with acos is right in case of small angles it is flat and acos give a fluctuations in angle by its precision lost when rounding.

But tan of a small angles is not flat, the direvative in 0 is 1 the problem of tan in a PI/2 (goes to infinity) but not for atan.

atan(close to infinity) -> PI/2

Did you notice the diffrence between
"atan" and "atan2" ?
atan2 - is a function of two arguments, the atan of the first argument divided by second. "atan2" function implemented in x86 FPU.

In case of atan2( x, y) let's play with this;

x = sin(angle);
y = cos(angle);
(i use "->" in term "tend" )
in a small angles
angle->0
the cos(angle) -> 1 and make fluctuations, but the sin(angle) -> angle

in a angles close to PI/2 or 3*PI/2 the
sin(angle) -> 1
and cos(angle) -> angle

So "atan2" use the precision of both sin and cos and don't introduce additional errors.

You can just test it :

float x = sin(angle);
float y = cos(angle);

float a_acos = acos(y);
float a_atan2 = atan2(x,y);

float acos_err = fabs(angle - a_acos);
float atan2_err = fabs(angle - a_atan2);

The "acos_err" always bigger or (not often) equal to "atan2_err". I got in most cases "atan2_err" - 0 but some time in angles that produce comparable sin and cos the error was near 1.0E-08 - 1.0E-12

By: minorlogic ( Michaele Norel )
RE: quaternion spherical interpolation SLERP  
2003-10-24 04:30

Hi Martin !
At the end i found the simpler formula of SLERP that not preserve the quaternion length.

look like
q =
q1*(tan(0.5*omega) - tan( (t - 0.5)*omega )) +
q2*(tan(0.5*omega) + tan( (t - 0.5)*omega ));

or

q =
(q1 + q2)*tan(0.5*omega) +
(q2 - q1)*tan( (t - 0.5)*omega );

And i tried to isolate the initialisation code from interpolation pass:
that what i get ( C++ code):
class SLERPer_tangent
{
xxquaternion q1, q2;
xxquaternion q_sum, q_dif;
double omega;
public:

SLERPer_tangent(const xxquaternion& q1_, const xxquaternion& q2_) :
q1( q1_ ), q2( q2_ )
{
q1.normalize();
q2.normalize();
float inner = inner_product(q1, q2);
if(inner < 0 ){
inner = -inner;
q2 = -q2;
}

if(inner > 0.99999f ) inner = 0.99999f;
omega = acos( inner );

q_sum = (q1 + q2)*(float)tan( omega*0.5 );
q_dif = (q2 - q1);
}

void Interpolate( float t, xxquaternion& qr)
{
/*float it = tan( (t - 0.5)*omega );
float s1 = tan_half_omega - it;
float s2 = tan_half_omega + it;
qr = q1*s1 + q2*s2 ;
rearange and simplify the equation
q1*s1 + q2*s2 = q1*(tan_half_omega - it) + q2*(tan_half_omega + it) =
(q1 + q2)*tan_half_omega + (q2 -q1)*it
*/
qr = q_sum + q_dif*tan( (t - 0.5)*omega );
}

};

But it perform the 1.3 -1.5 time faster than classic SLERP, i expected more ....

I think i stop for a time my reaserch in this , and will use the classic SLERP and LERP.