Advertisement

Inverse matrix doesn't actually inverse my rotation

Started by December 11, 2015 10:12 AM
3 comments, last by _WeirdCat_ 9 years, 2 months ago

so i try to yaw vector 0,0,1 by 90 degs with some matrix

then i try to inverse that matrix to get 0,0,1 again but i get 0,0,-1

i am not sure if code for matrix inverse is wrong or maybe the bug is here


	mat4 ahue;
	TEntity * p = new TEntity();
	p->rotation.yaw(cos(angle*imopi), sin(angle*imopi));
	p->rotation.DoRotation();
	p->rotation_mat.LoadGLMatrix(p->rotation.AIR_MATRIX);

vec3 zp = vec3(0.0, 0.0, 1.0);
vec3 zrot = p->rotation_mat* zp;

mat4 invrot;
invrot = p->rotation_mat;
invrot.Inverse();
vec3 irot = invrot * zrot;


mat4 res;
res = p->rotation_mat * invrot;

vec3 test = res * zp;

ALOG("BASE VEC: "+POINT_TO_TEXT(zp));
ALOG("ROTATED VEC: "+POINT_TO_TEXT(zrot));
ALOG("Inverted rotation VEC: "+POINT_TO_TEXT(irot));
ALOG("test VEC: "+POINT_TO_TEXT(test));

result:

12-11 11:06:09.779: V/WNB_LOG(11115): Called log but its disabled: BASE VEC: X 0 Y 0 Z 1

12-11 11:06:09.779: V/WNB_LOG(11115): Called log but its disabled: ROTATED VEC: X 1 Y 0 Z 6.12323e-17
12-11 11:06:09.779: V/WNB_LOG(11115): Called log but its disabled: Inverted rotation VEC: X 1.22465e-16 Y 0 Z -1
12-11 11:06:09.779: V/WNB_LOG(11115): Called log but its disabled: test VEC: X 1.22465e-16 Y 0 Z -1

where


template <class type> type  dp4(t4dpoint<type> vVector1, t4dpoint<type>  vVector2)
{
return ( (vVector1.x * vVector2.x) + (vVector1.y * vVector2.y) + (vVector1.z * vVector2.z) + (vVector1.w * vVector2.w) );
}

t3dpoint<T> operator*(const t3dpoint<T> p) const
{
t4dpoint<T> vertexPos;
vertexPos.x = p.x;
vertexPos.y = p.y;
vertexPos.z = p.z;
vertexPos.w = 1.0;

t3dpoint<T> vertexClip;
t4dpoint<T> matrow;

matrow.x = m[0]; matrow.y = m[1]; matrow.z = m[2]; matrow.w = m[3];
vertexClip.x = dp4(matrow, vertexPos);

matrow.x = m[4]; matrow.y = m[5]; matrow.z = m[6]; matrow.w = m[7];
vertexClip.y = dp4(matrow, vertexPos);

matrow.x = m[8]; matrow.y = m[9]; matrow.z = m[10]; matrow.w = m[11];
vertexClip.z = dp4(matrow, vertexPos);

return vertexClip;
}

and this is the the order i use: (i dont know which is that)


void Translate(T x, T y, T z)
{
Matrix44<T> tran(
 1, 0, 0, x,
 0, 1, 0, y,
 0, 0, 1, z,
 0, 0, 0, 1);
(*this) = (*this) * tran;
}

inverse code

[spoiler]



Matrix44<T>  Inverse()
{
 Matrix44<T> invOut;
  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 invOut;

    det = 1.0 / det;

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


		return invOut;

} 

[/spoiler]

A matrix times its inverse has to be the identity matrix. Your multiplication rot_mat * rot_mat_inv * (0,0,1) is not equal to (0,0,1) as it should, so either matrix * vector or matrix * matrix or your inversion code is broken. The best test for a matrix inversion routine is to multiply a matrix by its inverse and check that the result is equal to the identity matrix.

That aside, a matrix containing only a rotation is a so called orthogonal matrix, which means that its inverse is equal to its transpose (rows and columns swapped). So make sure the inversion of a rotation matrix results in the transposed rotation matrix.

I haven't checked your matrix inversion code, but another tidbit: if the determinant of a matrix is 0 it cannot be inverted. You're just returning some garbage right now, which is not a good way to handle such a (programmer) error...:)

Advertisement

invrot.Inverse();

That line of code computes the inverse and then throws it away.



invrot.Inverse();

That line of code computes the inverse and then throws it away.

Exactly this.

Also has as been pointed out, you should check if the determinant is 0. A better signature would be


bool Inverse(Matrix44 &invOut);

Then you can at least do:


mat4 invrot;

if(!p->rotation_mat.Inverse(invrot))
   // deal with being unable to invert

Interested in Fractals? Check out my App, Fractal Scout, free on the Google Play store.


invrot.Inverse();

That line of code computes the inverse and then throws it away.

LOL

x)

This topic is closed to new replies.

Advertisement