# Algorithm for 4x4 matrix inverse

This topic is 2121 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

## Recommended Posts

I need that to calculate the View matrix for OpenGL. I can do gaussian elimination on paper but I dont think I could implement that.

I found this though:

http://www.cg.info.hiroshima-cu.ac.jp/~miyazaki/knowledge/teche23.html

....is matrix inverse really such a complex thing to calculate? Isnt there a simpler algorithm for that? Thanks!

##### Share on other sites

I need that to calculate the View matrix for OpenGL. I can do gaussian elimination on paper but I dont think I could implement that.
I found this though:
http://www.cg.info.hiroshima-cu.ac.jp/~miyazaki/knowledge/teche23.html

You should probably look into LU-decomposition.

....is matrix inverse really such a complex thing to calculate? Isnt there a simpler algorithm for that? Thanks!

Yes, it is a complex thing.

##### Share on other sites

Since the view matrix is always orthonormal (or nearly so), can't you just use the transpose of the 3x3 part and do the translational part "by hand"? That's waaaaaaaay easier than doing an inverse.

Edited by samoth

##### Share on other sites
bool gluInvertMatrix(const double m[16], double invOut[16])
{
double inv[16], det;
int i;

inv[0] = m[5]  * m[10] * m[15] -
m[5]  * m[11] * m[14] -
m[9]  * m[6]  * m[15] +
m[9]  * m[7]  * m[14] +
m[13] * m[6]  * m[11] -
m[13] * m[7]  * m[10];

inv[4] = -m[4]  * m[10] * m[15] +
m[4]  * m[11] * m[14] +
m[8]  * m[6]  * m[15] -
m[8]  * m[7]  * m[14] -
m[12] * m[6]  * m[11] +
m[12] * m[7]  * m[10];

inv[8] = m[4]  * m[9] * m[15] -
m[4]  * m[11] * m[13] -
m[8]  * m[5] * m[15] +
m[8]  * m[7] * m[13] +
m[12] * m[5] * m[11] -
m[12] * m[7] * m[9];

inv[12] = -m[4]  * m[9] * m[14] +
m[4]  * m[10] * m[13] +
m[8]  * m[5] * m[14] -
m[8]  * m[6] * m[13] -
m[12] * m[5] * m[10] +
m[12] * m[6] * m[9];

inv[1] = -m[1]  * m[10] * m[15] +
m[1]  * m[11] * m[14] +
m[9]  * m[2] * m[15] -
m[9]  * m[3] * m[14] -
m[13] * m[2] * m[11] +
m[13] * m[3] * m[10];

inv[5] = m[0]  * m[10] * m[15] -
m[0]  * m[11] * m[14] -
m[8]  * m[2] * m[15] +
m[8]  * m[3] * m[14] +
m[12] * m[2] * m[11] -
m[12] * m[3] * m[10];

inv[9] = -m[0]  * m[9] * m[15] +
m[0]  * m[11] * m[13] +
m[8]  * m[1] * m[15] -
m[8]  * m[3] * m[13] -
m[12] * m[1] * m[11] +
m[12] * m[3] * m[9];

inv[13] = m[0]  * m[9] * m[14] -
m[0]  * m[10] * m[13] -
m[8]  * m[1] * m[14] +
m[8]  * m[2] * m[13] +
m[12] * m[1] * m[10] -
m[12] * m[2] * m[9];

inv[2] = m[1]  * m[6] * m[15] -
m[1]  * m[7] * m[14] -
m[5]  * m[2] * m[15] +
m[5]  * m[3] * m[14] +
m[13] * m[2] * m[7] -
m[13] * m[3] * m[6];

inv[6] = -m[0]  * m[6] * m[15] +
m[0]  * m[7] * m[14] +
m[4]  * m[2] * m[15] -
m[4]  * m[3] * m[14] -
m[12] * m[2] * m[7] +
m[12] * m[3] * m[6];

inv[10] = m[0]  * m[5] * m[15] -
m[0]  * m[7] * m[13] -
m[4]  * m[1] * m[15] +
m[4]  * m[3] * m[13] +
m[12] * m[1] * m[7] -
m[12] * m[3] * m[5];

inv[14] = -m[0]  * m[5] * m[14] +
m[0]  * m[6] * m[13] +
m[4]  * m[1] * m[14] -
m[4]  * m[2] * m[13] -
m[12] * m[1] * m[6] +
m[12] * m[2] * m[5];

inv[3] = -m[1] * m[6] * m[11] +
m[1] * m[7] * m[10] +
m[5] * m[2] * m[11] -
m[5] * m[3] * m[10] -
m[9] * m[2] * m[7] +
m[9] * m[3] * m[6];

inv[7] = m[0] * m[6] * m[11] -
m[0] * m[7] * m[10] -
m[4] * m[2] * m[11] +
m[4] * m[3] * m[10] +
m[8] * m[2] * m[7] -
m[8] * m[3] * m[6];

inv[11] = -m[0] * m[5] * m[11] +
m[0] * m[7] * m[9] +
m[4] * m[1] * m[11] -
m[4] * m[3] * m[9] -
m[8] * m[1] * m[7] +
m[8] * m[3] * m[5];

inv[15] = m[0] * m[5] * m[10] -
m[0] * m[6] * m[9] -
m[4] * m[1] * m[10] +
m[4] * m[2] * m[9] +
m[8] * m[1] * m[6] -
m[8] * m[2] * m[5];

det = m[0] * inv[0] + m[1] * inv[4] + m[2] * inv[8] + m[3] * inv[12];

if (det == 0)
return false;

det = 1.0 / det;

for (i = 0; i < 16; i++)
invOut[i] = inv[i] * det;

return true;
}


##### Share on other sites

To  explain samoth's post a bit:

The view matrix is the inverse of the camera's world transformation. This is usually a composition of translation and rotation, perhaps a scaling. Assuming column vector matrices, this looks like so:

C := T * R * S

A decomposition of C into a translation T, rotation R, and scaling S, as shown above, is relatively easy. It is even easier if scaling is known to not appear, because then decomposition is just extraction of values.

Then the inverse is

V := C-1S-1 * R-1 * T-1

because, due to associativity of the matrix product,
( T * R * ) * ( S-1 * R-1 * T-1) = T * ( R * ( S-1 ) * R-1 ) * T-1 = I

Looking at the properties of the particular matrices, the following correspondences can be used:
( S(sx, sy, sz) )-1 = S(1/sx, 1/sy, 1/sz)
( R )-1 = Rt
( T(tx, ty, tz) )-1 = T(-tx, -ty, -tz)

So the inversion of such a matrix C can be replaced by a decomposition and matrix products of usual translation, rotation, and scaling matrices.
Edited by haegarr

##### Share on other sites

Thank you for the answers! Thanks @haegarr for taking the time and explaining it. I implemented it and although I dont have the projection matrix yet, it seems OK.

• 21
• 11
• 9
• 17
• 13