Sign in to follow this  

4x4 Matrix Inversion

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

If you intended to correct an error in the post then please contact us.

Recommended Posts

Hi There all, I'm working on my first ever Graphics Engine and I'm stuck on the pretty hairy problem of Matrix inversion. I'm using the Building Blocks engine from the book Ultimate 3d Game Engine Design and Architecture as a base and then working off of that. The problem however, is that I have no Idea whatsoever what the code writer was doing when he wrote his matrix inversion code I really don't understand it. I've tried looking up the different methods but I still can't make any sense of it... Here's the code hopefully somebody will be able to tell me how it works or at least a basic idea of it. Thank you!

bool Matrix4x4::inverseMatrix(const Matrix4x4 &m)
{
   float tempMatrix[16] = {0};
   float d12, d13, d23, d24, d34, d41;

	d12 = matrix[2]  * matrix[7]  - matrix[3]  * matrix[6];
	d13 = matrix[2]  * matrix[11] - matrix[3]  * matrix[10];
	d23 = matrix[6]  * matrix[11] - matrix[7]  * matrix[10];
	d24 = matrix[6]  * matrix[15] - matrix[7]  * matrix[14];
	d34 = matrix[10] * matrix[15] - matrix[11] * matrix[14];
	d41 = matrix[14] * matrix[3]  - matrix[15] * matrix[2];

	tempMatrix[0] =   matrix[5] * d34 - matrix[9] * d24 + matrix[13] * d23;
	tempMatrix[1] = -(matrix[1] * d34 + matrix[9] * d41 + matrix[13] * d13);
	tempMatrix[2] =   matrix[1] * d24 + matrix[5] * d41 + matrix[13] * d12;
	tempMatrix[3] = -(matrix[1] * d23 - matrix[5] * d13 + matrix[9]  * d12);

	float determinant = matrix[0] * tempMatrix[0] + matrix[4] * tempMatrix[1] +
                       matrix[8] * tempMatrix[2] + matrix[12] * tempMatrix[3];

	if(determinant == 0.0)
	   {
		   Identity();
         return false;
	   }

	float invDeterminant = 1.0f / determinant;
	
	tempMatrix[0] *= invDeterminant;
	tempMatrix[1] *= invDeterminant;
	tempMatrix[2] *= invDeterminant;
	tempMatrix[3] *= invDeterminant;

	tempMatrix[4] = -(matrix[4] * d34 - matrix[8] * d24 + matrix[12] * d23) * invDeterminant;
	tempMatrix[5] =   matrix[0] * d34 + matrix[8] * d41 + matrix[12] * d13  * invDeterminant;
	tempMatrix[6] = -(matrix[0] * d24 + matrix[4] * d41 + matrix[12] * d12) * invDeterminant;
	tempMatrix[7] =   matrix[0] * d23 - matrix[4] * d13 + matrix[8]  * d12  * invDeterminant;

	d12 = matrix[0]  * matrix[5]  - matrix[1]  * matrix[12];
	d13 = matrix[0]  * matrix[9]  - matrix[1]  * matrix[8];
	d23 = matrix[4]  * matrix[9]  - matrix[5]  * matrix[8];
	d24 = matrix[4]  * matrix[13] - matrix[5]  * matrix[12];
	d34 = matrix[8]  * matrix[13] - matrix[9]  * matrix[12];
	d41 = matrix[12] * matrix[1]  - matrix[13] * matrix[0];

	tempMatrix[8]  =   matrix[7] * d34 - matrix[11] * d24 + matrix[15] * d23 * invDeterminant;
	tempMatrix[9]  = -(matrix[3] * d34 + matrix[11] * d41 + matrix[15] * d13) * invDeterminant;
	tempMatrix[10] =   matrix[3] * d24 + matrix[7]  * d41 + matrix[15] * d12 * invDeterminant;
	tempMatrix[11] = -(matrix[3] * d23 - matrix[7]  * d13 + matrix[11] * d12) * invDeterminant;
	tempMatrix[12] = -(matrix[6] * d34 - matrix[10] * d24 + matrix[14] * d23) * invDeterminant;
	tempMatrix[13] =   matrix[2] * d34 + matrix[10] * d41 + matrix[14] * d13 * invDeterminant;
	tempMatrix[14] = -(matrix[2] * d24 + matrix[6]  * d41 + matrix[14] * d12) * invDeterminant;
	tempMatrix[15] =   matrix[2] * d23 - matrix[6]  * d13 + matrix[10] * d12 * invDeterminant;

   matrix[0]  = tempMatrix[0];
   matrix[1]  = tempMatrix[1];
   matrix[2]  = tempMatrix[2];
   matrix[3]  = tempMatrix[3];

	matrix[4]  = tempMatrix[4];
	matrix[5]  = tempMatrix[5];
   matrix[6]  = tempMatrix[6];
   matrix[7]  = tempMatrix[7];

	matrix[8]  = tempMatrix[8];
	matrix[9]  = tempMatrix[9];
   matrix[10] = tempMatrix[10];
   matrix[11] = tempMatrix[11];

	matrix[12] = tempMatrix[12];
	matrix[13] = tempMatrix[13];
   matrix[14] = tempMatrix[14];
   matrix[15] = tempMatrix[15];

   return true;
}

Share this post


Link to post
Share on other sites
Quote:


bool Matrix4x4::inverseMatrix(const Matrix4x4 &m)
{
float tempMatrix[16] = {0};
float d12, d13, d23, d24, d34, d41;

d12 = matrix[2] * matrix[7] - matrix[3] * matrix[6];
d13 = matrix[2] * matrix[11] - matrix[3] * matrix[10];
d23 = matrix[6] * matrix[11] - matrix[7] * matrix[10];
d24 = matrix[6] * matrix[15] - matrix[7] * matrix[14];
d34 = matrix[10] * matrix[15] - matrix[11] * matrix[14];
d41 = matrix[14] * matrix[3] - matrix[15] * matrix[2];

tempMatrix[0] = matrix[5] * d34 - matrix[9] * d24 + matrix[13] * d23;
tempMatrix[1] = -(matrix[1] * d34 + matrix[9] * d41 + matrix[13] * d13);
tempMatrix[2] = matrix[1] * d24 + matrix[5] * d41 + matrix[13] * d12;
tempMatrix[3] = -(matrix[1] * d23 - matrix[5] * d13 + matrix[9] * d12);

float determinant = matrix[0] * tempMatrix[0] + matrix[4] * tempMatrix[1] +
matrix[8] * tempMatrix[2] + matrix[12] * tempMatrix[3];

if(determinant == 0.0)
{
Identity();
return false;
}



The above part calculates the determinant of the matrix. If a matrix has a determinant equal to zero the matrix is said to be Singular. Singular matrices are not invertible.

Quote:

float invDeterminant = 1.0f / determinant;

tempMatrix[0] *= invDeterminant;
tempMatrix[1] *= invDeterminant;
tempMatrix[2] *= invDeterminant;
tempMatrix[3] *= invDeterminant;

tempMatrix[4] = -(matrix[4] * d34 - matrix[8] * d24 + matrix[12] * d23) * invDeterminant;
tempMatrix[5] = matrix[0] * d34 + matrix[8] * d41 + matrix[12] * d13 * invDeterminant;
tempMatrix[6] = -(matrix[0] * d24 + matrix[4] * d41 + matrix[12] * d12) * invDeterminant;
tempMatrix[7] = matrix[0] * d23 - matrix[4] * d13 + matrix[8] * d12 * invDeterminant;

d12 = matrix[0] * matrix[5] - matrix[1] * matrix[12];
d13 = matrix[0] * matrix[9] - matrix[1] * matrix[8];
d23 = matrix[4] * matrix[9] - matrix[5] * matrix[8];
d24 = matrix[4] * matrix[13] - matrix[5] * matrix[12];
d34 = matrix[8] * matrix[13] - matrix[9] * matrix[12];
d41 = matrix[12] * matrix[1] - matrix[13] * matrix[0];

tempMatrix[8] = matrix[7] * d34 - matrix[11] * d24 + matrix[15] * d23 * invDeterminant;
tempMatrix[9] = -(matrix[3] * d34 + matrix[11] * d41 + matrix[15] * d13) * invDeterminant;
tempMatrix[10] = matrix[3] * d24 + matrix[7] * d41 + matrix[15] * d12 * invDeterminant;
tempMatrix[11] = -(matrix[3] * d23 - matrix[7] * d13 + matrix[11] * d12) * invDeterminant;
tempMatrix[12] = -(matrix[6] * d34 - matrix[10] * d24 + matrix[14] * d23) * invDeterminant;
tempMatrix[13] = matrix[2] * d34 + matrix[10] * d41 + matrix[14] * d13 * invDeterminant;
tempMatrix[14] = -(matrix[2] * d24 + matrix[6] * d41 + matrix[14] * d12) * invDeterminant;
tempMatrix[15] = matrix[2] * d23 - matrix[6] * d13 + matrix[10] * d12 * invDeterminant;

matrix[0] = tempMatrix[0];
matrix[1] = tempMatrix[1];
matrix[2] = tempMatrix[2];
matrix[3] = tempMatrix[3];

matrix[4] = tempMatrix[4];
matrix[5] = tempMatrix[5];
matrix[6] = tempMatrix[6];
matrix[7] = tempMatrix[7];

matrix[8] = tempMatrix[8];
matrix[9] = tempMatrix[9];
matrix[10] = tempMatrix[10];
matrix[11] = tempMatrix[11];

matrix[12] = tempMatrix[12];
matrix[13] = tempMatrix[13];
matrix[14] = tempMatrix[14];
matrix[15] = tempMatrix[15];

return true;
}



The last part does the actual inversion defined by the following equation:

A^-1 = (1 / det(A)) * adj(A)

Notice how we multiply the inverse determinant times the ajoint of matrix A. This definition of inversion goes along way to explain why if determinate is zero the matrix doesn't have an inverse.

The only other interesting thing about the above definition is the adjoint. You can Google or refer to a math text on how to calculate the ajoint of a matrix.


As a side note... I'm not entirely sure why you would want to invert a 4x4 matrix (presuming it is only to contain a 3D affine transformation). 4x4 transformation matrices are frequently non-invertible. For example consider the matrix with zero translation, that would produce a row with all zeros. It is known that a matrix with any row or column with all zeros has a determinant equal to zero, and is therefore non-invertible:

|1 0 0 0 |
|0 1 0 0 |
|0 0 1 1 |
|0 0 0 0 |

We would rather perform inversion by inverting the 3x3 linear transformation matrix in the upper left, and negating the translation portion in the bottom row of the 4x4 matrix.

Thats just my two-cents; I'm curious what other people have to say on that matter.

Share this post


Link to post
Share on other sites
Quote:
Original post by fpsgamer
As a side note... I'm not entirely sure why you would want to invert a 4x4 matrix (presuming it is only to contain a 3D affine transformation). 4x4 transformation matrices are frequently non-invertible. For example consider the matrix with zero translation, that would produce a row with all zeros. It is known that a matrix with any row or column with all zeros has a determinant equal to zero, and is therefore non-invertible:

|1 0 0 0 |
|0 1 0 0 |
|0 0 1 1 |
|0 0 0 0 |

We would rather perform inversion by inverting the 3x3 linear transformation matrix in the upper left, and negating the translation portion in the bottom row of the 4x4 matrix.

Thats just my two-cents; I'm curious what other people have to say on that matter.
This isn't quite right.

First of all, a 4x4 homogeneous transform matrix with a translation of (0, 0, 0) does not have a row of all zeros, but rather a row of [0 0 0 1] (remember that a 3D point represented in 4D homogeneous coordinates has a w value of 1).

Second of all, there are some cases in basic 3D math where a general 4x4 inversion is useful. Although you can special-case just about any inversion, a general inversion function can be useful for, say, inverting projection matrices for picking purposes (such inversions can also be special-cased, but IMO it's not really worth the trouble).

Also, although inversion of various affine transforms can certainly be special-cased, inverting the linear portion and negating the translation won't always give the correct results (in fact, it usually won't).

Share this post


Link to post
Share on other sites
I don't like the way this code tests for singularity. Testing (determinant == 0.0) is only asking for trouble. (fabsf(determinant) < EPSILON) is a far more stable alternative. I say so not only due to floating-point inaccuracy, but because this method of inversion will give stupid answers even for non-singular matrices, given a small enough determinant.

Share this post


Link to post
Share on other sites
re: why use a 4x4 matrice inversion in 3d space?

Identity matrices is 4x4

Quaternion Grassmann products are 4x4

in general, quaternion representation are either 2x2 complex number matrices, or 4x4 real number matricies

in my case - auto calibration of sensors (3D accelerometers) is a 6x6 for static, and 18x18 and greater for dynamic - not including movement and rotation tracking.


PS - after 2 or 3 years away from gamedev, I'm back :)

Share this post


Link to post
Share on other sites
Thanks for all your help! I was getting so confused whilst looking at the code and then looking at implementing a gaussian inversion method and then looking at the theory of cramer's rule. Hopefully I'll be able to implement my game engine in such a way that it actually runs. I'll keep the forums posted on my results.

-Jorick

Share this post


Link to post
Share on other sites

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

If you intended to correct an error in the post then please contact us.

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now

Sign in to follow this