Damn them quaternions
Author
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?
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?
What's the interpolation routine?
Author
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.
I won't be able to post the interpolation routine until tuesday unfortunately, but I'll update this post then.
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?
So is the problem you can't load the MD5 or you can't get your quaternion code to work right?
Author
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.
Author
The source code for my quaternions and the slerp routine (nicked from the net somewhere):
and the test program:
You'll notice that the slerped quaternion is not unit, which seems a bit odd to me, but that could be correct.
#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_HPPand 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.
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.
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.
Author
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.
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):
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.
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.
Author
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.
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.
This topic is closed to new replies.
Advertisement
Popular Topics
Advertisement
