🔒 Locked 🎮 Unity

Damn them quaternions

Started by rick_appleton Jul 31, 2005 at 4:20 AM 11 replies 3.3k views

Original Post

Jul 31, 2005 at 4:20 AM

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?

deffer

Jul 31, 2005 at 8:22 AM

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?

rick_appleton

Jul 31, 2005 at 8:44 AM

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.

Name_Unknown

Jul 31, 2005 at 7:39 PM

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?

rick_appleton

Aug 1, 2005 at 1:46 AM

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.

rick_appleton

Aug 2, 2005 at 2:11 AM

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.

deffer

Aug 2, 2005 at 5:17 AM

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.

rick_appleton

Aug 2, 2005 at 6:08 AM

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.

deffer

Aug 2, 2005 at 6:24 AM

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.

rick_appleton

Aug 2, 2005 at 5:19 PM

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.

deffer

Aug 2, 2005 at 7:29 PM

Actually, it is just perfect.
One note: rotating using interpolated quaternion is not the same as interpolating between two rotated vectors (using two base quaternions). There are many potential paths, that the interpolation could take. In the end, it takes only one of them, depending on method of interpolation used. But nowhere it is said that using slerp we are guaranteed to simulate (spherical)linear interpolation between two uninterpolated rotated-results!

It was kinda unintuitive to me at first, so I made a test program (using my particle engine, yay). I told the particles to flow from (0,0,0) in direction of the vector (1,0,0) rotated by a quaternion taken from interoplation between the two you gave, with random interpolation coefficient between 0 and 1.
In other "words":

const float r = randomf01();const quat q = slerp(m_q1, m_q2, r);velocity = vec(1,0,0) * q;

I took two screens to show the result. Hopefully, self-explanatory:
look first
look second

If you really want to have the results coplanar, you would have to explicitly use control points (they are like tangents in vector interpolation) and use ubic interpolation.
But I supose you don't care about it that much. The important thing is that it's calculated correctly.

bit64

Aug 2, 2005 at 7:45 PM

Rated Deffer 'extremely helpful' for this post.

rick_appleton

Aug 3, 2005 at 1:47 AM

So it appears the slerp is working correctly after all. I'd already rated you up before the last post, but I'd do it again if I could.

Thank you for taking the time to test that and to create those nice images. Truely appreciated.

So I guess the next step is testing my code with the doom 3 animation system again. See if it's still incorrect, or if the results are better now. I'll let you guys know the results.

Topic Locked

This topic has been locked by a moderator. New replies are not allowed.