Advertisement

Fastest way to inverse an orthogonal 4x4 matrix?

Started by September 24, 2004 04:28 AM
14 comments, last by snk_kid 20 years, 1 month ago
I know about the Kramer's rule (and have that working), the Gaussian method, and know that the inverse of a rotational 3x3 matrix is it's transpose. But I think I'm missing something. Since the matrix is orthogonal and I can calculate the 3x3 part of it by transposing it, shouldn't there be some easier method to calculate the remaining 3 numbers (the translate part)? I heard somewhere that they should be flipped from positive to negative, but that's not quite working. What am I missing? Here's an example... I have this matrix:
0  4  8  12
1  5  9  13
2  6  10 14
3  7  11 15

X:	0.819300	-0.196000	0.538800	-0.591540
Y:	0.301200	0.946800	-0.113700	-1.476260
Z:	-0.487900	0.255500	0.834700	-0.163410
W:	0.000000	0.000000	0.000000	1.000000
Using the Kramer's rule, the inverse of it is:
X:	0.819302	0.301249	-0.487825	0.849656
Y:	-0.195927	0.946701	0.255428	1.323417
Z:	0.538872	-0.113697	0.834705	0.287318
W:	0.000000	0.000000	0.000000	1.000000
But the center can be found by simply transposing the 3x3 part:
X:	0.819300	0.301200	-0.487900	-0.591540
Y:	-0.196000	0.946800	0.255500	-1.476260
Z:	0.538800	-0.113700	0.834700	-0.163410
W:	0.000000	0.000000	0.000000	1.000000
So the question is, how to calculate the m[12], m[13] and m[14] quickly? I suppose simply taking a part of the code from Kramer's rule for those 3 would do it, but that's cheap... Does anyone know a better way?
if you have 4x4 matrix A that acts as rotation matrix3x3 B with some translation V ,
that is, forward transform for X to X' is X'=A*X=B*X+V

So backward transform will be (matrix inverse is denoted by ~)
X = (~A)*X' = B'*X'+V' = (~B)*(X'-V) = (~B)*X'-(~B)*V
so
B' = (~B)
V' = - B'*V

that is, to compute inverse of 4X4 matrix with rotation and translation only(translation part is a last column) you need to:

transpose 3x3 matrix in top-left corner.
Multiply translation part by 3x3 matrix in top-left corner.
Negate translation part.
Advertisement
also, note that the lower row has three zeros and a one in it. if you know how to invert matrices, that should be enough to simplify things.
"3x3 part of it by transposing it"
/!\ Orthonormal (no scaling) is not orthogonal.

That's obvious if you use the linear properties.
Decompose your matrix by blocks as Mt*Mr. You first rotate by the 3x3 submatrix R then translate by T. Thus the inverse is
M^-1= Mr^-1 * Mt^-1
Consider R is made of columns X, Y, Z (the local axes).

 Mt             Mr            M 1  0  0  Tx    .  .  .  0    Xx Yx Zx Tx 0  1  0  Ty    .  R  .  0    Xy Yy Zy Ty 0  0  1  Tz    .  .  .  0    Xz Yz Zz Tz 0  0  0  1     0  0  0  1    0  0  0  1 Mt^-1          Mr^-1        M^-1 1  0  0 -Tx    .  .  .  0   Xx Xy Xz -X*Tx 0  1  0 -Ty    . R^-1.  0   Yx Yy Yz -Y*Ty 0  0  1 -Tz    .  .  .  0   Zx Zy Zz -Z*Tz 0  0  0  1     0  0  0  1   0  0  0  1


Comment :
In most typical cases, I would not even dare to inverse such a matrix. I would write a special routine named "transform vector by transposed" :

P' <- { X*(P-T), Y*(P-T), Z*(P-T) }

Edit : Dmytry and Eelco. Always the same club. You both beat me by a few seconds :) LOL !!!
"Coding math tricks in asm is more fun than Java"
Thanks! This works:

// Inverses an orthogonal matrixvoid CMatrix4x4::InvertOrthogonal(){	float fOut[16];	memcpy(fOut, fMat, sizeof(float) * 16);	fMat[1]  =  fOut[4];	fMat[2]  =  fOut[8];	fMat[4]  =  fOut[1];	fMat[6]  =  fOut[9];	fMat[8]  =  fOut[2];	fMat[9]  =  fOut[6];	fMat[12] = -(fOut[12] * fMat[0]  +  fOut[13] * fMat[4]  +  fOut[14] * fMat[8]);	fMat[13] = -(fOut[12] * fMat[1]  +  fOut[13] * fMat[5]  +  fOut[14] * fMat[9]);	fMat[14] = -(fOut[12] * fMat[2]  +  fOut[13] * fMat[6]  +  fOut[14] * fMat[10]);}
Quote: Original post by Charles B
"3x3 part of it by transposing it"
/!\ Orthonormal (no scaling) is not orthogonal.

Hehe. I'm a math noob, I can make such mistakes. x_X

At least it starts with the same "Ortho"! :D
Quote: That's obvious if you use the linear properties.
Decompose your matrix by blocks as Mt*Mr. You first rotate by the 3x3 submatrix R then translate by T. Thus the inverse is
M^-1= Mr^-1 * Mt^-1
Consider R is made of columns X, Y, Z (the local axes).

 Mt             Mr            M 1  0  0  Tx    .  .  .  0    Xx Yx Zx Tx 0  1  0  Ty    .  R  .  0    Xy Yy Zy Ty 0  0  1  Tz    .  .  .  0    Xz Yz Zz Tz 0  0  0  1     0  0  0  1    0  0  0  1 Mt^-1          Mr^-1        M^-1 1  0  0 -Tx    .  .  .  0   Xx Xy Xz -X*Tx 0  1  0 -Ty    . R^-1.  0   Yx Yy Yz -Y*Ty 0  0  1 -Tz    .  .  .  0   Zx Zy Zz -Z*Tz 0  0  0  1     0  0  0  1   0  0  0  1


Comment :
In most typical cases, I would not even dare to inverse such a matrix. I would write a special routine named "transform vector by transposed" :

P' <- { X*(P-T), Y*(P-T), Z*(P-T) }

Edit : Dmytry and Eelco. Always the same club. You both beat me by a few seconds :) LOL !!!

You guys are all awesome. ;) Much appreciated~
Advertisement
Quote: Original post by MichaelMook
Thanks! This works:


OK but pls, optimize it a bit now that it works.

Use temporary floats instead of copying the full matrix on the stack. Also beware memcpy is sometimes very slow, in debug for instance. Even John Carmack fell in the trap once. He wondered why his code was so slow and it took him several hours to see it came from the std C library implementation. It's somthing I knew before he wrote his first books, I have always checked what the compilers do. I am a speedy code addict, low level biased :)

You need 9 temporary variables. With some luck, since the FPU has 8 registers, only one such variable will be on the stack.

Write something like :

float f4=M[4], f8=M[8], etc ...;

M[1]=f4;
"Coding math tricks in asm is more fun than Java"
Obviously for transformation matrices the other presentations here, especially Charles B's, are THE appropriate way to do this if you absolutely must have the inverse.

But, suppose you have a general, non-transformation matrix...

Just a note. For matrices around 4x4 and above, Cramer's rule is going to be the slowest method! It's an O((n+1)!)) operation. For a general matrix 3x3 or larger, use Gaussian Elimination with partial pivoting or one of the more advanced direct method. Very large matrices (maybe 1000x1000 or larger) require an iterative method, since roundoff will kill Gaussian Elimination.

[edited by grhodes_at_work to correct the operation count for Cramer's Rule. Its far, far worse than O(n3).]

[Edited by - grhodes_at_work on September 25, 2004 7:32:48 PM]
Graham Rhodes Moderator, Math & Physics forum @ gamedev.net
Quote: Original post by Charles B
Quote: Original post by MichaelMook
Thanks! This works:


OK but pls, optimize it a bit now that it works.

Use temporary floats instead of copying the full matrix on the stack. Also beware memcpy is sometimes very slow, in debug for instance. Even John Carmack fell in the trap once. He wondered why his code was so slow and it took him several hours to see it came from the std C library implementation. It's somthing I knew before he wrote his first books, I have always checked what the compilers do. I am a speedy code addict, low level biased :)

You need 9 temporary variables. With some luck, since the FPU has 8 registers, only one such variable will be on the stack.

Write something like :

float f4=M[4], f8=M[8], etc ...;

M[1]=f4;

Thanks, I'll keep that in mind. I didn't notice any change in framerate, but I will keep it in mind for the future when I have more happening in the scene and will be able to judge more precisely.

Quote: Just a note. For matrices around 4x4 and above, Cramer's rule is going to be the slowest method!

Yup, I saw your other posts ;). I'm not really using it though. I also ended up making two other versions of the above function, as it seems I was getting eventual rounding errors by not using the determinant. And the other one merely flips the translation coordinates (seems to work nicely for billboarding of the sun and requires the least number of calculations).
Sorry i haven't read the whole thread but i've been using Cramer's rule expressed differently to get the inverse matrix explicitly for 4x4 matrix and below:

M-1 = (1 / det(M)) * adj(M)

where det = determinate &
adj = adjoint.

i thought this was the fastest way for matrices upto 4x4? anything bigger than 4x4 use Gaussian elimination.

and for orthogonal matrices:

M-1 = MT

[Edited by - snk_kid on September 24, 2004 5:40:23 PM]

This topic is closed to new replies.

Advertisement