# Matrix to Quaternion calculation instable

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

## Recommended Posts

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.

##### Share on other sites

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

##### Share on other sites

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;
}
}


1. 1
2. 2
Rutin
21
3. 3
A4L
15
4. 4
5. 5

• 13
• 26
• 10
• 11
• 44
• ### Forum Statistics

• Total Topics
633740
• Total Posts
3013621
×