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.