Maths - Conversion Matrix to Quaternion - Forum Discussion

By: Nobody/Anonymous - nobody
RE: Matrix to Quaternion error.  
2006-01-14 15:08

Hi Martin, 
 
I just stumbled over your matrix-quat conversion page, which is very good imho. 
 
Here are a couple of suggestions: 
 
1) Because quaternions cannot represent a reflection component, the matrix must be special orthogonal. For a special orthogonal matrix, 1 + trace is always positive. The case switch is not needed, just do sqrt() on the 4 trace variants and you are done, ie: 
 
quaternion.w = sqrt( max( 0, 1 + m00 + m11 + m22 ) ) / 2; 
quaternion.x = sqrt( max( 0, 1 + m00 - m11 - m22 ) ) / 2; 
quaternion.y = sqrt( max( 0, 1 - m00 + m11 - m22 ) ) / 2; 
quaternion.z = sqrt( max( 0, 1 - m00 - m11 + m22 ) ) / 2; 
 
The max( 0, ... ) is just a safeguard against rounding error. 
 
 
2) You can make the conversion deal with positive scale as well if you relax the assumption that det(matrix)=1. For this, the norm of the quaternion must be known a priori (the cube root of the volume spanned by the matrix axes). 
 
absQ2 = det( matrix )^(1/3) 
 
quaternion.w = sqrt( max( 0, absQ2 + m00 + m11 + m22 ) ) / 2;  
 
...etc 
 
cheers 
 

By: Martin Baker - martinbaker
RE: Matrix to Quaternion error.  
2006-01-15 08:44

"For a special orthogonal matrix, 1 + trace is always positive"
 
Thanks very much Ill update the comments on the web page. 
 
I guess the advantage of the case switch is that we only need to do a single sqrt instead of 4 sqrts which might be a performance issue. 
 
The other issue that interests me is sensitivity to rounding errors. In other words, the matrix contains redundant information so rounding errors may de-orthogonalise the matrix, if this happens which method is most likely to cancel out the errors? 
 
I dont know how best to tackle this? say, for instance, rounding errors introduce an error "delta e" to one element of the matrix, but we dont know which element has the error. Which method is most likely to cancel out this error? Which method is most likely to produce a normalised quaternion? 
 
Thanks, 
 
Martin

By: Nobody/Anonymous - nobody
RE: Matrix to Quaternion error.  
2006-01-15 09:33

 
Hi martin, 
 
Sorry, I forgot to paste the part that recovers the signs as well. You need to do this after the sqrt(). 
 
Q.x = _copysign( Q.x, m21 - m12 ) 
Q.y = _copysign( Q.y, m02 - m20 ) 
Q.z = _copysign( Q.z, m20 - m02 ) 
 
Depdending on your convention, the signs may be reverse. 
 
My Intel datasheet says, a sqrt() is 24 .. 40 cycles, depending on the current prescision setting in the float-control register. A division is comparable to a sqrt() in terms of cycles.  
 
So you win the branchlessness which is better IMHO (a mispredicted branch can cost you another 50 cycles, and since we're branching on random data, it will mispredict 50% of time, the worst case). 
 
Another pitfall you might watch out is whether _copysign() is implemented as intrinsic and not as a library call. 
 
However, _copysign is really just the transfer of the highest bit of the floating point number, so you could write a macro that does that. 
 
 
The most robust code will have an expression that consider all redundant values at once. For instance, if you know the squared length of the quaternion must be equal to the length any row- or column-vectors, then instead of picking out a signle one, better calculate something that involves all values at the same time. 
 
 
Christian 

By: Nobody/Anonymous - nobody
RE: Matrix to Quaternion error.  
2006-02-06 16:21

Hello Martin, 
 
I think I have a couple of corrections to your webpage on the Matrix to quaternion convertions. 
First, I like Christians method because it is  
more straight-forward and you aviod all the IF statements. But I think he made a typo in the last equation: 
 
I think it should be: 
 
quaternion.w = sqrt( max( 0, 1 + m00 + m11 + m22 ) ) / 2; 
quaternion.x = sqrt( max( 0, 1 + m00 - m11 - m22 ) ) / 2; 
quaternion.y = sqrt( max( 0, 1 - m00 + m11 - m22 ) ) / 2; 
quaternion.z = sqrt( max( 0, 1 - m00 - m11 + m22 ) ) / 2; 
Q.x = _copysign( Q.x, m21 - m12 ) 
Q.y = _copysign( Q.y, m02 - m20 ) 
Q.z = _copysign( Q.z, m10 - m01 ) <--correction 
 
 
On the other method with the IF statements, I think the signs should be reversed on some of the qw's. It's not a big deal. All it means is that the sign on qw is reversed when qw**2 < 0.25*epsilon. Here is how I think it should be: 
 
IF (T > epsilon) 
 
S = 0.5 / sqrt(T) 
 
W = 0.25 / S 
 
X = ( m21 - m12 ) * S 
 
Y = ( m02 - m20 ) * S 
 
Z = ( m10 - m01 ) * S 
 
ELSE 
 
if (m00 > m11)&(m00 > m22)) { 
S = sqrt( 1.0 + m00 - m11 - m22 ) * 2; 
qx = 0.25 * S; 
qy = (m01 + m10 ) / S; 
qz = (m02 + m20 ) / S; 
qw = (m21 - m12 ) / S; <--correction 
 
} else if (m11 > m22)) { 
S = sqrt( 1.0 + m11 - m00 - m22 ) * 2; 
qx = (m01 + m10 ) / S; 
qy = 0.25 * S; 
qz = (m12 + m21 ) / S; 
qw = (m02 - m20 ) / S; 
 
} else { 
S = sqrt( 1.0 + m22 - m00 - m11 ) * 2; 
qx = (m02 + m20 ) / S; 
qy = (m12 + m21 ) / S; 
qz = 0.25 * S; 
qw = (m10 - m01 ) / S; <--correction 
} 
 
 
Thanks for a good webpage. 
-Stig 
 

By: Martin Baker - martinbaker
RE: Matrix to Quaternion error.  
2006-02-07 09:30

Thanks very much, I have made these corrections to the page. 
 
Martin