There is a potential issue with the use of 'copysign' as described here.
By: Nobody/Anonymous  nobody
RE: Matrix to Quaternion error.
20060114
15:08 
Hi Martin,
I just stumbled over
your matrixquat 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.
20060115
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
deorthogonalise 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.
20060115
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 floatcontrol 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
columnvectors, 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.
20060206
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 straightforward 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

This site may have errors. Don't use for critical systems.