Jump to content

  • Log In with Google      Sign In   
  • Create Account

Matrix to Quaternion calculation instable


Old topic!
Guest, the last post of this topic is over 60 days old and at this point you may not reply in this topic. If you wish to continue this conversation start a new topic.

  • You cannot reply to this topic
2 replies to this topic

#1 RPTD   Members   -  Reputation: 316

Like
0Likes
Like

Posted 13 July 2014 - 06:55 AM

For animation purpose I move between matrix and quaterions in different places and for this I use the tracing method found on the internet:

const double trace = a11 + a22 + a33 + 1.0;
if( trace > 0.0001 ) {
    const double s = 0.5 / sqrt( trace );
    return decQuaternion( ( a32 - a23 ) * s, ( a13 - a31 ) * s, ( a21 - a12 ) * s, 0.25 / s );
}else if( a11 > a22 && a11 > a33 ){
    const double s = 2.0 * sqrt( 1.0 + a11 - a22 - a33 );
    return decQuaternion( 0.25 * s, ( a12 + a21 ) / s, ( a13 + a31 ) / s, ( a23 - a32 ) / s );
}else if( a22 > a33 ){
    const double s = 2.0 * sqrt( 1.0 + a22 - a11 - a33 );
    return decQuaternion( ( a12 + a21 ) / s, 0.25 * s, ( a23 + a32 ) / s, ( a13 - a31 ) / s );
}else{
    const double s = 2.0 * sqrt( 1.0 + a33 - a11 - a22 );
    return decQuaternion( ( a13 + a31 ) / s, ( a23 + a32 ) / s, 0.25 * s, ( a12 - a21 ) / s );
}

The matrix is in row major order and quaterions are in the (x,y,z,w) format.

 

If I do for example a small sweep (from [0,175°,0] to [0,185°,0]) across the XZ plane (hence with Y axis fixed to [0,1,0] where I'm using DX coordinate system) around the backwards pole (0,180°,0) I end up with a slight twiching of the camera near the [0,180°,0] point. I tracked it down to the calculated quaterion to be slightly off near the point where you use any other than the first if-case. Augmenting the step value to 0.0001 did help in some cases but not in this one here. I even went all the way up to 0.01 in which case the slight twiching just moved a bit further away from the problem point.

 

I also do not think the quaterion-to-matrix is the culprit since this code there does not use any if-cases and thus should be stable. Furthermore tweaking the above mentioned value does modify the error behavior so it has to be this code causing troubles. I can at the time being cheat around the problem but I'm looking for a proper solution.

 

So my question is what other possibility is there to calculate a quaterion from a matrix which is stable? Is there a known problem with the trace based method used here that doesn't work around the backwards point? I'm concerned more about an error free solution than the fastest solution on earth.


Life's like a Hydra... cut off one problem just to have two more popping out.
Leader and Coder: Project Epsylon | Drag[en]gine Game Engine

Sponsor:

#2 Aressera   Members   -  Reputation: 1484

Like
0Likes
Like

Posted 13 July 2014 - 02:25 PM

Try this conversion method, it is the same one used in the OGRE engine. Note: totally untested in this form.

/// Create a new quaternion from a 3x3 orthonormal rotation matrix.
/// Quaternion = {a, b, c, d}, Matrix3D = column major
Quaternion( const Matrix3D<T>& m )
{
	// Algorithm in Ken Shoemake's article in 1987 SIGGRAPH course notes
	// article "Quaternion Calculus and Fast Animation".
	T fTrace = m.x.x + m.y.y + m.z.z;
	
	if ( fTrace > T(0) )
	{
		// |w| > 1/2, may as well choose w > 1/2
		T fRoot = math::sqrt(fTrace + 1.0f);  // 2w
		a = 0.5f*fRoot;
		fRoot = 0.5f/fRoot;  // 1/(4w)
		
		b = (m.y.z - m.z.y)*fRoot;
		c = (m.z.x - m.x.z)*fRoot;
		d = (m.x.y - m.y.x)*fRoot;
	}
	else
	{
		// |w| <= 1/2
		Index nextIndex[3] = { 1, 2, 0 };
		Index i = 0;
		
		if ( m.y.y > m.x.x )
			i = 1;
		
		if ( m.z.z > m[i][i] )
			i = 2;
		
		Index j = nextIndex[i];
		Index k = nextIndex[j];
		
		T fRoot = math::sqrt( m[i][i] - m[j][j] - m[k][k] + T(1) );
		
		T* apkQuat[3] = { &b, &c, &d };
		*apkQuat[i] = T(0.5)*fRoot;
		
		fRoot = T(0.5)/fRoot;
		
		a = (m[j][k] - m[k][j])*fRoot;
		*apkQuat[j] = (m[i][j] + m[j][i])*fRoot;
		*apkQuat[k] = (m[i][k] + m[k][i])*fRoot;
	}
}

Edited by Aressera, 13 July 2014 - 02:26 PM.


#3 RPTD   Members   -  Reputation: 316

Like
0Likes
Like

Posted 26 July 2014 - 06:28 PM

Sorry for not replying sooner but the forum failed to send me a notification about you having posted so I assumed nobody answered until I decided to still check back in case the forum soft is broken (as it is).

 

So to your post... I could be wrong but is this not the same as my version just written in a non-un-rolled version? Maybe I'm missing something but it looks similar to me.

 

 

Try this conversion method, it is the same one used in the OGRE engine. Note: totally untested in this form.

/// Create a new quaternion from a 3x3 orthonormal rotation matrix.
/// Quaternion = {a, b, c, d}, Matrix3D = column major
Quaternion( const Matrix3D<T>& m )
{
	// Algorithm in Ken Shoemake's article in 1987 SIGGRAPH course notes
	// article "Quaternion Calculus and Fast Animation".
	T fTrace = m.x.x + m.y.y + m.z.z;
	
	if ( fTrace > T(0) )
	{
		// |w| > 1/2, may as well choose w > 1/2
		T fRoot = math::sqrt(fTrace + 1.0f);  // 2w
		a = 0.5f*fRoot;
		fRoot = 0.5f/fRoot;  // 1/(4w)
		
		b = (m.y.z - m.z.y)*fRoot;
		c = (m.z.x - m.x.z)*fRoot;
		d = (m.x.y - m.y.x)*fRoot;
	}
	else
	{
		// |w| <= 1/2
		Index nextIndex[3] = { 1, 2, 0 };
		Index i = 0;
		
		if ( m.y.y > m.x.x )
			i = 1;
		
		if ( m.z.z > m[i][i] )
			i = 2;
		
		Index j = nextIndex[i];
		Index k = nextIndex[j];
		
		T fRoot = math::sqrt( m[i][i] - m[j][j] - m[k][k] + T(1) );
		
		T* apkQuat[3] = { &b, &c, &d };
		*apkQuat[i] = T(0.5)*fRoot;
		
		fRoot = T(0.5)/fRoot;
		
		a = (m[j][k] - m[k][j])*fRoot;
		*apkQuat[j] = (m[i][j] + m[j][i])*fRoot;
		*apkQuat[k] = (m[i][k] + m[k][i])*fRoot;
	}
}

Life's like a Hydra... cut off one problem just to have two more popping out.
Leader and Coder: Project Epsylon | Drag[en]gine Game Engine




Old topic!
Guest, the last post of this topic is over 60 days old and at this point you may not reply in this topic. If you wish to continue this conversation start a new topic.



PARTNERS