Decomposition of 3D (4x4) Transformation Matrices
Kirk Parsons
I am having a hard time getting the math to work, in a case
where I am trying to decompose a transformation matrix into
a combination of scale, rotation and translation matrices.
I've got code that works most of the time, but seems to be
incorrectly handling some special case, or singularity.
To anyone who can help, thanks in advance.
- Kirk Parsons
An example of a test case that causes my code to fail is as
follows. I start out with the VRML 2.0 transforms:
T1 = translation -0.15 -2.49 9.503
R1 = rotation 0.2164 -0.00945 -0.9763 -3.056
S = scale -1 -1 -1.08
R2 = scaleOrientation -0.1827 0.5362 0.8241 -0.7841
And the following matrix was built by applying the transforms
in the order (which is, I believe the correct interpretation
of the VRML above):
M = T1 R1 R2 S -R2
Using the convention of the transformation matrix pre-multiplying
a column vector, this gives us the matrix:
0.902868 0.0875462 0.420965 -0.15
-0.0793841 0.996177 -0.0369111 -2.49
0.456399 8.03066e-05 -0.979028 9.503
0 0 0 1
When I decompose the matrix, I get the following decomposition
(which I believe has an incorrect rotation):
translation -0.15 -2.49 9.503
rotation -.194848 .308740 .930975 2.96179
scale -1.01478 -1.00002 -1.06633
When I rebuild the matrix off this data, I should get the same matrix
that I started with (within roundoff error). Instead I get the
following matrix:
0.92197 0.2858 0.324798 -.15
-0.0478429 0.794785 -0.645062 -2.49
0.421213 -0.535381 -0.784199 9.503
0 0 0 1
I am fairly confident that I am computing this matrix correctly, and
that my error lies in the way that I decompose the matrix into
translation, rotation and scale transformations.
I have laid out some pseudo code for two functions, DecomposeMatrix
and GetRotationFromMatrix. Of these two functions, it seems more likely
that the error is in DecomposeMatrix. DecomposeMatrix takes a
transformation
matrix as input, and produces a translation, rotation and scale
transformations. The resulting transformations would be applied to
a point in local coordinates (r), to get the point in world coordinates
(r')
in the following equation:
r' = T*R*S*r
where T is the translation matrix, R is the rotation matrix and S is the
scale matrix.
If anyone out there has written similar code, I have a couple of
specific
questions that would be helpful to get answered.
1) It would be useful to run the matrix above through working algorithm,
to
see what the correct result is. Anyone have the right answer for the
matrix
above?
2) In the algorithm specified below, the part that I least understand
(and
therefore most likely to contain the error), is where I look for the
right-handedness of the incoming matrix. All the test cases that cause
my
code to fail so far have involved this piece of code. I believe that it
is also true that my code fails when the original matrix is built with
one
or more negative scale factors.
If I don't check for right handedness, then the GetRotationFromMatrix
fails,
but on the other hand, I'm not at all sure that I am handling this case
correctly. Has anyone out there identified what to do for left handed
matrices? Or for transform matrices that include negative scales?
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void DecomposeMatrix(AgMatrix matrix)
{
.
.
.
//
---------------------------------------------------------------------------
// Find magnitude of scale factors, by computing sqrt's of column
// vectors, then set translation vector:
//
// Note: rotation matrix has elements:
// r00 r01 r02 r03
// r10 r11 r12 r13
// r20 r21 r22 r23
// 0 0 0 1
rotation_matrix = matrix;
scale->x = (float) sqrt(r00*r00 + r10*r10 + r20*r20);
scale->y = (float) sqrt(r01*r01 + r11*r11 + r21*r21);
scale->z = (float) sqrt(r02*r02 + r12*r12 + r22*r22);
translate->x = r03;
translate->y = r13;
translate->z = r23;
//
---------------------------------------------------------------------------
// Factor out scale and remove translation to make rotation_matrix a
true
// rotation matrix.
col_1 = first column vector of rot matrix;
col_2 = second column vector of rot matrix;
col_3 = third column vector of rot matrix;
col_1 = col_1/scale->x;
col_2 = col_2/scale->y;
col_3 = col_3/scale->z;
r03 = 0;
r13 = 0;
r23 = 0;
//
---------------------------------------------------------------------------
// Check to see if rotation matrix is right handed.
//
// If dot prod if of cross_1x2 and col_3 is negative, then the matrix
is left
// handed. If this is the case, we negate both the scale and the
rotation
// matrix, which should leave us with the equivalent rotation
CrossProduct(&cross_1x2, &col1, &col2);
if ( cross_1x2.x*col_3.x + cross_1x2.y*col_3.y + cross_1x2.z*col_3.z <
0)
{
scale->x = -scale->x;
scale->y = -scale->y;
scale->z = -scale->z;
rotation_matrix = - rotation_matrix;
}
GetRotationFromMatrix(&rotation_matrix, &rotation_vector, &rotation);
}
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
The following pseudo code is a description of the algorithm in my
GetRotationFromMatrix function, which calculates a rotation axis, plus
rotation
angle for a given rotation matrix. It seems more likely that my error is
in the
DecomposeMatrix function above, than in this function, but it is
included here
for reference.
This algorithm is from Graphics Gems I, "Rotation Tools", with the
addition of
some special case code for a rotation angle near 0. In the particular
case
where our test data fails, the special case code is not being executed.
void AgGetRotationFromMatrix(rotation_matrix, rotation_vector, rotation)
{
cos_theta = 0.5 * (r00 + r11 + r22 - 1);
// arccos(0.999961923) == 0.5 degrees
double fSpecialCaseCk = 0.999961923;
if (cosTheta <= -fSpecialCaseCk) rotation = PI;
else if (cosTheta >= fSpecialCaseCk) rotation = 0.0f;
else rotation = acos(cos_theta);
if (rotation == PI) {
diag1 = 0.5 * (r00 + 1);
diag2 = 0.5 * (r11 + 1);
diag3 = 0.5 * (r22 + 1);
if ((diag1>=diag2) && (diag1>=diag3)) {
// diag1 is largest diagonal element
rotation_vector->x = sqrt(diag1);
rotation_vector->y = 0.5 * r01/rotation_vector->x;
rotation_vector->z = 0.5 * r02/rotation_vector->x;
} else if ((diag2>=diag1) && (diag2>=diag3)) {
// diag2 is largest diagonal element
rotation_vector->y = (float)sqrt(diag2);
rotation_vector->x = 0.5 * r10/rotation_vector->y;
rotation_vector->z = 0.5 * r12/rotation_vector->y;
} else {
// diag3 is largest diagonal element
rotation_vector->z = (float)sqrt(diag3);
rotation_vector->x = 0.5 * r20/rotation_vector->z;
rotation_vector->y = 0.5 * r21/rotation_vector->z;
}
} else if (rotation == 0) {
rotation_vector->x = 1.0f;
rotation_vector->y = 0.0f;
rotation_vector->z = 0.0f;
} else {
// NOTE: in the test case data listed above, this is the branch
// of the if statement that is executed.
rotation_vector->x = r21 - r12;
rotation_vector->y = r02 - r20;
rotation_vector->z = r10 - r01;
}
Normalize(rotation_vector);
return eE_OK;
}
=======================================================
Kirk Parsons
Avatar Software Designer
Attic Graphics, Inc.
ki...@mindspring.com
David Eberly
In article <346CF1E7...@mindspring.com>,
Kirk Parsons <ki...@mindspring.com> wrote:
>I am having a hard time getting the math to work, in a case
>where I am trying to decompose a transformation matrix into
>a combination of scale, rotation and translation matrices.
I used my singular value decomposition code on your matrix
of interest and had no problems with the decomposition. The
scale values were (1.080183,1.000023,1.000010).
I think your DecomposeMatrix is not correct. It is not
enough to unitize the columns of a matrix to get an orthonormal
matrix. The columns need to also be mutually orthogonal (the
dot products of pairs of columns must be zero).
Dave Eberly
ebe...@cs.unc.edu