Hi, I''m looking for a good algo to calculate the inverse of a 4x4 matrix. I didn''t find anything good with google. I mean, I can have a lot of information about calculating the inverse of a matrix on paper, but that''s mostly impossible to implement as C code So how do you guys calculate matrix inverses in your engines/games/projects? Any code or information would be useful. Thanks, Flo

code snipeds:

//#define _float double#ifndef _float#define _float	float#endifclass MATRIX{// Last row is 0, 0, 0, 1, else change class functions( remove // and /* ), frustum changes last row	union{		_float m[4*4]; // m[x][y]= [((x-1)*4) + (y-1)]		struct {	// mxy !!!(in definition/here reversed)=>			_float	m11, m12, m13, m14,	//	m11 m21 m31 m41									m21, m22, m23, m24,	//	m12 m22 m32 m42					m31, m32, m33, m34,	//	m13 m23 m33 m43					m41, m42, m43, m44;	//	m14 m24 m34 m44		};	};...

_float MATRIX::GetM3Determinant(const MATRIX3 &matrix){	return (	matrix.x.x * (matrix.y.y * matrix.z.z - matrix.y.z * matrix.z.y) -				matrix.y.x * (matrix.x.y * matrix.z.z - matrix.x.z * matrix.z.y) +				matrix.z.x * (matrix.x.y * matrix.y.z - matrix.x.z * matrix.y.y)	);}MATRIX::MATRIX3 MATRIX::GetSubMatrix(long x, long y){	MATRIX3 tmp;	long xoffset= 0;	for(long i=0; i < 4; i++)		if(i != x){			long yoffset= 0;			for(long j= 0; j < 4; j++)				if(j != y){					*( ((_float *) &tmp) + xoffset*3 + yoffset)= m[i*4 + j];					yoffset++;				}			xoffset++;		}	return tmp;}_float MATRIX::GetDeterminant(){	_float result= 0;	long i= 1;	MATRIX3 tmp;	for(long n= 0; n < 4; n++, i*=-1){		tmp= GetSubMatrix(0, n);		result+= m[n] * GetM3Determinant(tmp) * i;	}	return result;}MATRIX MATRIX::GetCorrectInverse(){	MATRIX inverse;	MATRIX3 tmp;	_float m4determinant= GetDeterminant();	if(fabs(m4determinant) < bias)		inverse= MATRIX(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0);// Matrix has no inverse	else{		for( long i=0; i < 4; i++)			for(long j=0; j < 4; j++){				long sign= 1-((i+j)%2)*2;					tmp= GetSubMatrix(i, j);					inverse.m[i + j*4]= GetM3Determinant(tmp) * sign / m4determinant;			}	}	return inverse;}

i hope you can read the important part, but i think i will edit my post in few minutes...

edit:

class VECTOR3D{public:	_float x, y, z;...   	struct MATRIX3{		VECTOR3D x,y,z;	};

uhm i dont think it will compile, i have not posted overloaded operators etc...


edit2:

_float MATRIX::GetDeterminant(){//of a 4x4 Matrix!	_float result= 0;	long i= 1;	MATRIX3 tmp;	for(long n= 0; n < 4; n++, i*=-1){		tmp= GetSubMatrix(0, n);		result+= m[n] * GetM3Determinant(tmp) * i;	}	return result;}

edit3: added source tags (this board eats them...


T2k

[edited by - T2k on September 13, 2003 6:36:55 AM]

[edited by - T2k on September 13, 2003 6:40:18 AM]

[edited by - T2k on September 13, 2003 6:43:13 AM]

That''s a bit of code
I think I''ll have to read this many times to understand.

So do you obtain submatrices until you get a 2x2 matrix and calculate the determinants like this, or what?

|a b|
|c d| det = ab - cd

I''ll try to implement this in my matrix class then

quote:

---

Original post by ZMaster
|a b|
|c d| det = ab - cd

---

det= ad - cb


T2k

Please, help me remember what inversed matrix are used for ? Inverse kynematic maybe ? and what else ?

Just an observation. While it is important to have a good general method for finding an inverse, always keep in mind that there are several special cases that can greatly simplify things.

For example, if a transformation matrix is built only from rotations, then the inverse of the matrix is guaranteed to be the transpose.

You can also often get an inverse matrix by inverting parameters. For example, a matrix given by
glTranslatef(10.0,10.0,10.0) is guaranteed to be the inverse of
glTranslatef(-10.0, -10.0, -10.0).


[edited by - UnIcron on September 15, 2003 12:20:12 PM]

MV, in OpenGL is the most common use of the inverse if you want to set the camera. Set it as an object but invert before sendinf it to OpenGL.

In that case, isn''t it easier to store the camera''s position and orientation as a vector and a quaternion ?

---

SaM3d!, a cross-platform API for 3d based on SDL and OpenGL.
The trouble is that things never get better, they just stay the same, only more so. -- (Terry Pratchett, Eric)

Well, thanks. I got a good routine now. I''m building submatrices (first 3x3, then 2x2) and just calculate the inverse of those (like T2k suggested).

Thank for the tips UnIcron. Should I implement some different routines for performance reason? I''ll think about that - maybe I should test it

Basicly: The matrix multiplied by it''s inverse equals the identity matrix.
I don''t have a practical example in mind now, but you must sometimes use a matrix inverse when converting coordinates to different spaces. I''m coding some math classes for my engine, and I wanted to implement an Invert routine.

@rodzilla: That''s easier, sure, but it always depends on what camera you want to set up. I finished my 3rd person camera class some minutes ago. I use matrices in it very frequently. There is a function Perspective(), that calculates a perspective projection matrix and uses it in OpenGL (just like gluPerspective) for example. Moreover I have a matrix that takes the camera transforms (rotation, translation). I''s multiplied with the current modelview matrix every frame, so I could easily write a function, that automatically tracks an object (like glLookAt, but more automatic) just by making some changes to the camera transform matrix.

float Determinant4f(const float m[16])
{
return
m[12]*m[9]*m[6]*m[3]-
m[8]*m[13]*m[6]*m[3]-
m[12]*m[5]*m[10]*m[3]+
m[4]*m[13]*m[10]*m[3]+
m[8]*m[5]*m[14]*m[3]-
m[4]*m[9]*m[14]*m[3]-
m[12]*m[9]*m[2]*m[7]+
m[8]*m[13]*m[2]*m[7]+
m[12]*m[1]*m[10]*m[7]-
m[0]*m[13]*m[10]*m[7]-
m[8]*m[1]*m[14]*m[7]+
m[0]*m[9]*m[14]*m[7]+
m[12]*m[5]*m[2]*m[11]-
m[4]*m[13]*m[2]*m[11]-
m[12]*m[1]*m[6]*m[11]+
m[0]*m[13]*m[6]*m[11]+
m[4]*m[1]*m[14]*m[11]-
m[0]*m[5]*m[14]*m[11]-
m[8]*m[5]*m[2]*m[15]+
m[4]*m[9]*m[2]*m[15]+
m[8]*m[1]*m[6]*m[15]-
m[0]*m[9]*m[6]*m[15]-
m[4]*m[1]*m[10]*m[15]+
m[0]*m[5]*m[10]*m[15];
}

BOOL GenerateInverseMatrix4f(float i[16], const float m[16])
{
float x=Determinant4f(m);
if (x==0) return FALSE;

i[0]= (-m[13]*m[10]*m[7] +m[9]*m[14]*m[7] +m[13]*m[6]*m[11]
-m[5]*m[14]*m[11] -m[9]*m[6]*m[15] +m[5]*m[10]*m[15])/x;
i[4]= ( m[12]*m[10]*m[7] -m[8]*m[14]*m[7] -m[12]*m[6]*m[11]
+m[4]*m[14]*m[11] +m[8]*m[6]*m[15] -m[4]*m[10]*m[15])/x;
i[8]= (-m[12]*m[9]* m[7] +m[8]*m[13]*m[7] +m[12]*m[5]*m[11]
-m[4]*m[13]*m[11] -m[8]*m[5]*m[15] +m[4]*m[9]* m[15])/x;
i[12]=( m[12]*m[9]* m[6] -m[8]*m[13]*m[6] -m[12]*m[5]*m[10]
+m[4]*m[13]*m[10] +m[8]*m[5]*m[14] -m[4]*m[9]* m[14])/x;
i[1]= ( m[13]*m[10]*m[3] -m[9]*m[14]*m[3] -m[13]*m[2]*m[11]
+m[1]*m[14]*m[11] +m[9]*m[2]*m[15] -m[1]*m[10]*m[15])/x;
i[5]= (-m[12]*m[10]*m[3] +m[8]*m[14]*m[3] +m[12]*m[2]*m[11]
-m[0]*m[14]*m[11] -m[8]*m[2]*m[15] +m[0]*m[10]*m[15])/x;
i[9]= ( m[12]*m[9]* m[3] -m[8]*m[13]*m[3] -m[12]*m[1]*m[11]
+m[0]*m[13]*m[11] +m[8]*m[1]*m[15] -m[0]*m[9]* m[15])/x;
i[13]=(-m[12]*m[9]* m[2] +m[8]*m[13]*m[2] +m[12]*m[1]*m[10]
-m[0]*m[13]*m[10] -m[8]*m[1]*m[14] +m[0]*m[9]* m[14])/x;
i[2]= (-m[13]*m[6]* m[3] +m[5]*m[14]*m[3] +m[13]*m[2]*m[7]
-m[1]*m[14]*m[7] -m[5]*m[2]*m[15] +m[1]*m[6]* m[15])/x;
i[6]= ( m[12]*m[6]* m[3] -m[4]*m[14]*m[3] -m[12]*m[2]*m[7]
+m[0]*m[14]*m[7] +m[4]*m[2]*m[15] -m[0]*m[6]* m[15])/x;
i[10]=(-m[12]*m[5]* m[3] +m[4]*m[13]*m[3] +m[12]*m[1]*m[7]
-m[0]*m[13]*m[7] -m[4]*m[1]*m[15] +m[0]*m[5]* m[15])/x;
i[14]=( m[12]*m[5]* m[2] -m[4]*m[13]*m[2] -m[12]*m[1]*m[6]
+m[0]*m[13]*m[6] +m[4]*m[1]*m[14] -m[0]*m[5]* m[14])/x;
i[3]= ( m[9]* m[6]* m[3] -m[5]*m[10]*m[3] -m[9]* m[2]*m[7]
+m[1]*m[10]*m[7] +m[5]*m[2]*m[11] -m[1]*m[6]* m[11])/x;
i[7]= (-m[8]* m[6]* m[3] +m[4]*m[10]*m[3] +m[8]* m[2]*m[7]
-m[0]*m[10]*m[7] -m[4]*m[2]*m[11] +m[0]*m[6]* m[11])/x;
i[11]=( m[8]* m[5]* m[3] -m[4]*m[9]* m[3] -m[8]* m[1]*m[7]
+m[0]*m[9]* m[7] +m[4]*m[1]*m[11] -m[0]*m[5]* m[11])/x;
i[15]=(-m[8]* m[5]* m[2] +m[4]*m[9]* m[2] +m[8]* m[1]*m[6]
-m[0]*m[9]* m[6] -m[4]*m[1]*m[10] +m[0]*m[5]* m[10])/x;

return TRUE;
}

[edited by - circuit on September 16, 2003 2:08:36 PM]