note: when we started this discussion the code on the page was:
If the trace of the matrix is greater than zero, then the result is:
S = 0.5 / sqrt(T)
W = 0.25 / S
X = ( m21 - m12 ) * S
Y = ( m02 - m20 ) * S
Z = ( m10 - m01 ) * S
If the trace of the matrix is less than or equal to zero then identify which major diagonal element has the greatest value.
if ((m00 > m11)&(m00 > m22)) {
S = sqrt( 1.0 + m00 - m11 - m22 ) * 2; // S=4*qx
qw = (m21 - m12) / S;
qx = 0.25 * S;
qy = (m01 + m10) / S;
qz = (m02 + m20) / S;
} else if (m11 > m22) {
S = sqrt( 1.0 + m11 - m00 - m22 ) * 2; // S=4*qy
qw = (m02 - m20) / S;
qx = (m01 + m10) / S;
qy = 0.25 * S;
qz = (m12 + m21) / S;
} else {
S = sqrt( 1.0 + m22 - m00 - m11 ) * 2; // S=4*qz
qw = (m10 - m01) / S;
qx = (m02 + m20) / S;
qy = (m12 + m21) / S;
qz = 0.25 * S;
}
| RE: Matrix to Quaternion error. By: Ethan Tira-Thompson (ejtttje) - 2009-06-17 18:09 |
|
Hi Martin, I was looking at the matrixToQuaternion page, and I think the "Better code" and "C++ code" sections may have a problem with numerical stability as Tr + 1 may approach 0 in some cases. The test for greater-than-epsilon helps, but a better solution may be found in the implementation from the Wild Magic package by Dave Eberly: http://www.geometrictools.com/LibFoundation/Mathematics/Wm4Quaternion.inl (See "FromRotationMatrix" function) There, they compare Tr > 0 (instead of Tr+1 > 0) for the initial conditional, thus avoiding getting close to 0 in the 0.5 / sqrt(Tr + 1) which follows it. Do you have any thoughts, does this make sense to you? The implications of falling into the 'else' clause are the main concern here, but that seems safe, and appears to work in practice. Thanks! |
|
Hi Ethan, Thanks, yes I agree with what you say, I guess there are two issues: 1) danger of division by zero or square root of a negative number. if m11 and m22 are very close to m00 it may pass the (m00 > m11)&(m00 > m22) test but the square root may still give zero or negative due to floating point errors. I don't know how likely this is but its better to keep clear. 2) accuracy of dividing by (and square rooting) a very small number. with floating point numbers, dividing small numbers by small numbers should be reasonably accurate but at the extreme it would loose accuracy. So given these two issues its seems to me the most accurate method would be to calculate all 4 major diagonals and then choose the one that is greatest rather than just being above '0'. The Dave Eberly code (using a Ken Shoemake algorithm according to the comment) only switches between two cases. I think I need to make some changes to make absolutely sure the code is resilient but I'm not quite sure how much complexity its worth including to maximise the accuracy. I'll have to give this some more thought. Martin |
RE: Matrix to Quaternion error.
By: Ethan Tira-Thompson (ejtttje) - 2009-06-19 01:24 |
|
> 1) danger of division by zero or square root of a negative number. > 2) accuracy of dividing by (and square rooting) a very small number. I think because we pick the biggest diagonal element in the 'else', and subtract the others, we can never get a negative trace in the 'else' cases. Assuming the following proof holds, then because of the +1 in the sqrt calls, the divisors won't become "small", nor negative (more precisely, won't become less than 1) A quick attempt at a proof that the chosen sqrt argument (1 + X ± Y ± Z) is always >= 1: A) Assume the biggest trace (X) is positive a) Assume second biggest (Y) is positive 1) if third biggest (Z) is positive, the trace has to be positive (X + Y + Z case) 2) if trace with third is negative, know X - Y > 0, so X - Y - Z will also be > 0 because Z negative b) Assume second biggest (Y) is negative 1) third biggest must also be negative, then X - Y - Z > 0 because X > 0 and we're subtracting negatives B) Assume biggest (X) is negative a) second and third must be more negative: X - Y - Z will have to be > 0 (actually, will even be >= |X|) So I don't think there's a downside in avoiding the <0 case in the initial positive-trace condition since the 'else' conditions will also avoid getting close to 0... |
|
Yes, I do like the approach of using the largest diagonal like this: float tr1 = 1.0 + m00 - m11 - m22
float tr2 = 1.0 - m00 + m11 - m22
float tr3 = 1.0 - m00 - m11 + m22
if (tr1 > tr2)&(tr1 > tr3)) {
float S = sqrt(tr1) * 2; // S=4*qx
qw = (m21 - m12) / S;
qx = 0.25 * S;
qy = (m01 + m10) / S;
qz = (m02 + m20) / S;
} else if (tr2 > tr1)&(tr2 > tr3)) {
float S = sqrt(tr2) * 2; // S=4*qy
qw = (m02 - m20) / S;
qx = (m01 + m10) / S;
qy = 0.25 * S;
qz = (m12 + m21) / S;
} else {
float S = sqrt(tr3) * 2; // S=4*qz
qw = (m10 - m01) / S;
qx = (m02 + m20) / S;
qy = (m12 + m21) / S;
qz = 0.25 * S;
}
I don't know if this is valid but I wonder if this would also be more
resilient to a de-orthogonalised matrix. For example, if before doing
this conversion, we have some calculations and small floating point
errors have made the matrix slightly de-orthogonalised. It seems to me
that the largest elements of the matrix would be the most accurate? and
if we choose the largest denominator then we must have the largest
numerator. |
