Jump to content
  • Advertisement
Sign in to follow this  
rick_appleton

Damn them quaternions

This topic is 4759 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

If you intended to correct an error in the post then please contact us.

Recommended Posts

I've recently created a Doom3 animation/mesh loader and although the animation is going well with my quaternions, I can't get them to blend correctly (currently I'm jumping from one frame to the next). I've searched on the net for some source code, but it seems (s)lerping quaternions doesn't exactly produce the results I expect. The situation: - A point at (1,0,0) - One quaternion which rotates a point -90 degrees around the y-axis (resulting in above point moving to (0,0,1) - Another quaternion which rotates a point +90 degrees around the z-axis (resulting in the starting point moving to (0,1,0)) - A slerped quaternion slerping the above two quaternions with factor 0.5 (for now) As I thought of quaternion this would result in a rotation of 90 degrees around the (0,-1,1) axis (normalized) thus resulting in the original point moving to (0, sqrt(2.f)/2.f, sqrt(2.f)/2.f). But the code I've seen results in something like (x, sqrt(2.f)/2.f, sqrt(2.f)/2.f) I haven't checked the exact result, but x>0. Is this correct?

Share this post


Link to post
Share on other sites
Advertisement
Rotation is not supposed to change the length of a vector being rotated, and here it definitely did. Unless x is a result of rounding errors (but if it was, you wouldn't have brought that up).

What's the interpolation routine?

Share this post


Link to post
Share on other sites
Thank you for the reply. This at least confirms that I'm doing something wrong.

I won't be able to post the interpolation routine until tuesday unfortunately, but I'll update this post then.

Share this post


Link to post
Share on other sites
Well if doesn't work with lerp, than you did something wrong loading it. Lerp will work, I have done it. Slerp will work too, I also have done it and it works fine with the MD5.

So is the problem you can't load the MD5 or you can't get your quaternion code to work right?

Share this post


Link to post
Share on other sites
The loading works fine, and animating the models based on the info in the anim file also seems to work correctly, however I am not currently blending between different frames of animation, simply switching from one to the next. So I suspect my quaternion code to be at fault. This is also confirmed in the test case described above.

Share this post


Link to post
Share on other sites
The source code for my quaternions and the slerp routine (nicked from the net somewhere):

#ifndef INC_MATH_QUATERNION_HPP
#define INC_MATH_QUATERNION_HPP

#include "Math/matrix.hpp"

template <class T = float>
class Quaternion
{
public:

//Quaternion
// -default constructor
// -creates a new quaternion with all parts equal to zero
Quaternion(void)
{
x = y = z = w = 0;
}


//Quaternion
// -constructor
// -parametes : x, y, z, w elements of the quaternion
// -creates a new quaternion based on the elements passed in
Quaternion(T wi, T xi, T yi, T zi)
{
w = wi;
x = xi;
y = yi;
z = zi;
}


//Quaternion
// -constructor
// -parameters : vector/array of four elements
// -creates a new quaternion based on the elements passed in
Quaternion(T v[4])
{
w = v[0];
x = v[1];
y = v[2];
z = v[3];
}


//Quaternion
// -copy constructor
// -parameters : const quaternion q
// -creates a new quaternion based on the quaternion passed in
Quaternion(const Quaternion<T>& q)
{
w = q.w;
x = q.x;
y = q.y;
z = q.z;
}

//Quaternion
// -constructor
// -parameters : vector + float for axis and angle
// -creates a new quaternion based on the elements passed in
Quaternion( float fAngle, vec3f axis )
{
w = cosf( fAngle/2);
x = axis.x * sinf( fAngle/2 );
y = axis.y * sinf( fAngle/2 );
z = axis.z * sinf( fAngle/2 );
}

//~Quaternion
// -destructor
// -deleted dynamically allocated memory
~Quaternion()
{
}


//operator=
// -parameters : q1 - Quaternion object
// -return value : Quaternion
// -when called on quaternion q2 sets q2 to be an object of q3
Quaternion<T> operator = (const Quaternion<T>& q)
{
w = q.w;
x = q.x;
y = q.y;
z = q.z;

return (*this);
}

//operator+
// -parameters : q1 - Quaternion object
// -return value : Quaternion
// -when called on quaternion q2 adds q1 + q2 and returns the sum in a new quaternion
Quaternion<T> operator + (const Quaternion<T>& q)
{
return Quaternion(w+q.w, x+q.x, y+q.y, z+q.z);
}

//operator-
// -parameters : q1- Quaternion object
// -return values : Quaternion
// -when called on q1 subtracts q1 - q2 and returns the difference as a new quaternion
Quaternion<T> operator - (const Quaternion<T>& q)
{
return Quaternion(w-q.w, x-q.x, y-q.y, z-q.z);
}


//operator*
// -parameters : q1 - Quaternion object
// -return values : Quaternion
// -when called on a quaternion q2, multiplies q2 *q1 and returns the product in a new quaternion
Quaternion<T> operator * (const Quaternion<T>& q) const
{
return Quaternion(
w*q.w - x*q.x - y*q.y - z*q.z,
w*q.x + x*q.w + y*q.z - z*q.y,
w*q.y + y*q.w + z*q.x - x*q.z,
w*q.z + z*q.w + x*q.y - y*q.x );
}

//operator*
// -parameters : f - float
// -return values : Quaternion
Quaternion<T> operator * ( float f ) const
{
return Quaternion( w*f, x*f, y*f, z*f );
}

//operator/
// -parameters : q1 and q2- Quaternion objects
// -return values : Quaternion
// -divide q1 by q2 and returns the quotient q1
Quaternion<T> operator / (Quaternion<T>& q)
{
return ((*this) * (q.inverse()));
}


//operator+=
// -parameters : q1- Quaternion object
// -return values : Quaternion
// -when called on quaternion q3, adds q1 and q3 and returns the sum as q3
Quaternion<T>& operator += (const Quaternion<T>& q)
{
w += q.w;
x += q.x;
y += q.y;
z += q.z;

return (*this);
}


//operator-=
// -parameters : q1- Quaternion object
// -return values : Quaternion
// -when called on quaternion q3, subtracts q1 from q3 and returns the difference as q3
Quaternion<T>& operator -= (const Quaternion<T>& q)
{
w -= q.w;
x -= q.x;
y -= q.y;
z -= q.z;

return (*this);
}


//operator*=
// -parameters : q1- Quaternion object
// -return values : Quaternion
// -when called on quaternion q3, multiplies q3 by q1 and returns the product as q3
Quaternion<T>& operator *= (const Quaternion<T>& q)
{
T w_val = w*q.w - x*q.x - y*q.y - z*q.z;
T x_val = w*q.x + x*q.w + y*q.z - z*q.y;
T y_val = w*q.y + y*q.w + z*q.x - x*q.z;
T z_val = w*q.z + z*q.w + x*q.y - y*q.x;

w = w_val;
x = x_val;
y = y_val;
z = z_val;

return (*this);
}


//operator/=
// -parameters : q1- Quaternion object
// -return values : quaternion
// -when called on quaternion q3, divides q3 by q1 and returns the quotient as q3
Quaternion<T>& operator /= (Quaternion<T>& q)
{
(*this) = (*this)*q.Inverse();
return (*this);
}


//operator!=
// -parameters : q1 and q2- Quaternion objects
// -return value : bool
// -determines if q1 and q2 are not equal
bool operator != (const Quaternion<T>& q)
{
return (w!=q.w || x!=q.x || y!=q.y || z!=q.z) ? true : false;
}

//operator==
// -parameters : q1 and q2- Quaternion objects
// -return value : bool
// -determines if q1 and q2 are equal
bool operator == (const Quaternion<T>& q)
{
return (w==q.w && x==q.x && y==q.y && z==q.z) ? true : false;
}

//norm
// -parameters : none
// -return value : T
// -when called on a quaternion object q, returns the norm of q
T Norm()
{
return (w*w + x*x + y*y + z*z);
}

//magnitude
// -parameters : none
// -return value : T
// -when called on a quaternion object q, returns the magnitude q
T Magnitude()
{
return sqrt(Norm());
}

//scale
// -parameters : s- a value to scale q1 by
// -return value: quaternion
// -returns the original quaternion with each part, w,x,y,z, multiplied by some scalar s
Quaternion<T> Scale(T s)
{
return Quaternion(w*s, x*s, y*s, z*s);
}

// -parameters : none
// -return value : quaternion
// -when called on a quaternion object q, returns the inverse of q
Quaternion<T> Inverse()
{
return Conjugate().Scale(1/Norm());
}

//conjugate
// -parameters : none
// -return value : quaternion
// -when called on a quaternion object q, returns the conjugate of q
Quaternion<T> Conjugate()
{
return Quaternion(w, -x, -y, -z);
}

//UnitQuaternion
// -parameters : none
// -return value : quaternion
// -when called on quaterion q, takes q and returns the unit quaternion of q
Quaternion<T> MakeUnit()
{
return (*this).Scale(1/(*this).Magnitude());
}

// -parameters : vector of type T
// -return value : void
// -when given a 3D vector, v, rotates v by this quaternion

vec3f QuatRotation(vec3f v)
{
Quaternion <T> qv(0, v.x, v.y, v.z);
Quaternion <T> qm = (*this) * qv * (*this).Inverse();

return vec3f(qm.x, qm.y, qm.z);
}

matrix44 MakeMatrix( void ) const
{
matrix44 out;
out.mMatrix[0] = 1-2*y*y-2*z*z;
out.mMatrix[4] = 2*x*y-2*w*z;
out.mMatrix[8] = 2*x*z+2*w*y;
out.mMatrix[12] = 0;

out.mMatrix[1] = 2*x*y+2*w*z;
out.mMatrix[5] = 1-2*x*x-2*z*z;
out.mMatrix[9] = 2*y*z-2*w*z;
out.mMatrix[13] = 0;

out.mMatrix[2] = 2*x*z-2*w*y;
out.mMatrix[6] = 2*y*z-2*w*z;
out.mMatrix[10] = 1-2*x*x-2*y*y;
out.mMatrix[14] = 0;

out.mMatrix[3] = 0;
out.mMatrix[7] = 0;
out.mMatrix[11] = 0;
out.mMatrix[15] = 1;
return out;
}

float w,x,y,z;

friend Quaternion<T> QuaternionSlerp(const Quaternion<T>& q1, const Quaternion<T>& q2, const float t )
{
float to1[4];
float cosom = q1.x * q2.x + q1.y * q2.y + q1.z * q2.z + q1.w * q2.w;
if (cosom<0)
{
cosom = -cosom;
to1[0] = -q2.x;
to1[1] = -q2.y;
to1[2] = -q2.z;
to1[3] = -q2.w;
}
else
{
to1[0] = q2.x;
to1[1] = q2.y;
to1[2] = q2.z;
to1[3] = q2.w;
}

// Slerp
/*float omega = acos(cosom);
float sinom = sin(omega);
float scale0 = sin((1.f - t)*omega)/sinom;
float scale1 = sin(t*omega)/sinom;*/


// Lerp
float scale0 = (1.f - t);
float scale1 = t;


Quaternion<float> qr;
qr.x = scale0*q1.x + scale1*to1[0];
qr.y = scale0*q1.y + scale1*to1[1];
qr.z = scale0*q1.z + scale1*to1[2];
qr.w = scale0*q1.w + scale1*to1[3];

return qr;
}
/*{
float s = q1.x * q2.x + q1.y * q2.y + q1.z * q2.z + q1.w * q2.w;
if (s>1)
s=1;
else if (s<-1)
s=-1;

float phi = acos(s);

Quaternion<float> qr;

float sin_phi = sinf(phi);
float sin_phit = sinf(phi*t);
float sin_phi_phit = sinf(phi-phi*t);

qr = q1*(sin_phi_phit/sin_phi) + q2*(sin_phit/sin_phi);
return qr;
}/*
{
float sq1, sq2;

float cosom = q1.x * q2.x + q1.y * q2.y + q1.z * q2.z + q1.w * q2.w;

Quaternion<T> q2b, qr;

if( cosom < 0.0f)
{
cosom = -cosom;
q2b.x = -q2.x;
q2b.y = -q2.y;
q2b.z = -q2.z;
q2b.w = -q2.w;
}
else
{
q2b = q2;
}

if( (1.0f + cosom) > 1E-5)
{
if( (1.0f - cosom) > 1E-5)
{
float om = acos(cosom);
float rsinom = (1.0f / sin(om));

sq1 = sinf( (1.0f - t) * om) * rsinom;
sq2 = sinf(t * om) * rsinom;
}
else
{
sq1 = (1.0f - t);
sq2 = t;
}

qr.w = sq1 * q1.w + sq2 * q2b.w;
qr.x = sq1 * q1.x + sq2 * q2b.x;
qr.y = sq1 * q1.y + sq2 * q2b.y;
qr.z = sq1 * q1.z + sq2 * q2b.z;
}
else
{
const float PI = 3.14159265358979323846f;

sq1 = sinf( (1.0f - t) * 0.5f * PI);
sq2 = sinf(t * 0.5f * PI);

qr.w = sq1 * q1.w + sq2 * q1.z;
qr.x = sq1 * q1.x + sq2 * q1.y;
qr.y = sq1 * q1.y + sq2 * q1.x;
qr.z = sq1 * q1.z + sq2 * q1.w;
}

return qr;
}*/

};
#endif // INC_MATH_QUATERNION_HPP



and the test program:

vec3f start(1,0,0), outA, outB, outC, slerped;
Quaternion<float> q1(-90*3.14/180, vec3f(0,1,0));
Quaternion<float> q2(90*3.14/180, vec3f(0,0,1));
outA = q1.QuatRotation(start); // (0,0,1)
outB = q2.QuatRotation(start); // (0,1,0)
vec3f tmp(0,-1,1);
tmp.Normalize();
Quaternion<float> q3(90*3.14/180, tmp);
outC = q3.QuatRotation(start); // (0, sqrt(2)/2, sqrt(2)/2);

Quaternion<float> q3 = QuaternionSlerp(q1,q2, 0.5);
q3.MakeUnit();
slerped = q3.QuatRotation(start); // (0.3..., 0.65, 0.65);



You'll notice that the slerped quaternion is not unit, which seems a bit odd to me, but that could be correct.

Share this post


Link to post
Share on other sites
Well, how to say this...
Did you noticed that the actual slerp code is commented out?
That there is lerp used instead?

Guess it's time to find another implementation.
I created mine based on the right way. The whole site is worth reading imo.

Share this post


Link to post
Share on other sites
Yes, I did :) Should have added that above. There are three implementations, and neither of them seems to produce the expected result. I'll take a look at that site.

Share this post


Link to post
Share on other sites
Ok, I don't know what's wrong with the commented version.

Here's my code for reference (I use it for a long time so i'm pretty sure it works):

// user expects result angle betw.(0, pi/2).
INL float angleQ1(const xquat & q1, const xquat & q2)
{
return 2.0f * asinf( (float) halve(lenq( subq(q1,q2) )) );
};
// user expects result angle betw.(pi/2, pi).
INL float angleQ2(const xquat & q1, const xquat & q2)
{
return 2.0f * asinf( (float) halve(lenq( addq(q1,q2) )) );
};

// angle between two unit-length quaternions.
// it is more stable alias for: acos(dot_q(q1, q2)).
INL float angleQ(const xquat & q1, const xquat & q2)
{
if ( isneg(dot_q(q1, q2)) ) {
// second quarter.
return F_PI - angleQ2(q1, q2);
} else {
// first quarter.
return angleQ1(q1, q2);
};
};

// fAng is a hint value equal to: angleQ(q1, q2).
// fSinAngOverAng is a hint value equal to: fSinOverX(fAng).
INL xquat slerp(const xquat & q1, const xquat & q2, float t, float fAng, float fSinAngOverAng)
{
const float t1 = 1.0f - t;
const float w1 = t1 * fSinOverX(t1 * fAng) / fSinAngOverAng;
const float w2 = t * fSinOverX(t * fAng) / fSinAngOverAng;

return add_weightq(q1, q2, w1, w2);
};
// shorten version.
INL xquat slerp(const xquat & q1, const xquat & q2, float t)
{
const float fAng = angleQ(q1, q2);
return slerp(q1, q2, t, fAng, fSinOverX(fAng));
};



addq, subq, dot_q, isneg, add_weightq - should be self-explanatory, I think.
halve - divides by 2.

The way is taken completely from that site.

Share this post


Link to post
Share on other sites
Thanks for that deffer. I've tried your version, and at least the slerp returns a unit quaternion now, but with the example I gave above, I'm getting the same resulting point (ie x!=0). Could you check that the result I'm expecting is actually what I should be getting? If so, then one of my other member functions is at fault, but at least I know I should keep looking. Unfortunately I'm a bit pressed for time the next two weeks, so I don't know when I'll have the time.

Just as a summary:
- Quaternion1 rotates -90 degrees around the y-axis
- Quaternion2 rotates +90 degrees around the z-axis
- Quaternion3 is a slerp of Quaternion1 and Quaternion2 with factor 0.5
- Rotate (1,0,0) over Quaternion3
- The result is ???

If this is undefined due to angle being pi/2, then please use:
Quaternion2 is +90 degrees rotation around axis (0,-1,1) with unit length.

I've tried this second one as well, and although the result seems better, x!=0.

Share this post


Link to post
Share on other sites
Sign in to follow this  

  • Advertisement
×

Important Information

By using GameDev.net, you agree to our community Guidelines, Terms of Use, and Privacy Policy.

We are the game development community.

Whether you are an indie, hobbyist, AAA developer, or just trying to learn, GameDev.net is the place for you to learn, share, and connect with the games industry. Learn more About Us or sign up!

Sign me up!