1

I am trying to get the inverse of a 4x4 square matrix for opengl, so it is a column major matrix. Also I would like to avoid discussions about making my code into loops, it's quite challenging enough to follow without loops.

I wrote some code that can do most basic matrix operations but matrix inverse was kinda rough, I found this answer: inverting a 4x4 matrix

I implemented it in my code

void mat4inverse(mat4* out, const mat4* in)
{
    double inv[16], det;
    int i;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    if (det == 0)
        printf("oh oh!");

    det = 1.0 / det;
    //printf("\n\ndeterminant:%f\n\n",det);

    for (i = 0; i < 16; i++)
        out->m[i] = inv[i] * det;
}

I create a matrix that looks like this:

[   1.000000    3.000000    4.000000    9.000000    ]
[   5.000000    6.000000    7.000000    2.000000    ]
[   8.000000    8.000000    8.000000    9.000000    ]
[   0.000000    0.000000    0.000000    1.000000    ]

after passing it to my transpose function that I listed above the results are:

[   -1.000000   1.000000    -0.375000   10.375000   ]
[   2.000000    -3.000000   1.625000    -26.625000  ]
[   -1.000000   2.000000    -1.125000   15.125000   ]
[   0.000000    0.000000    0.000000    1.000000    ]

I am not sure but that seems a bit wrong to me. This is for opengl rendering and my matrix is a 1D 16 double array.

Am I getting correct results?

additional info Here is my math over on code review: https://codereview.stackexchange.com/questions/98692/4d-matrix-math-library-for-use-with-opengl

Community
  • 1
  • 1
user1610950
  • 1,837
  • 5
  • 33
  • 49
  • 1
    Your inverted matrix looks fine to me. Try it out here http://www.bluebit.gr/matrix-calculator/calculate.aspx – Banex Jul 31 '15 at 22:29
  • 1
    both code and the result look fine – Sung Jul 31 '15 at 23:29
  • thank you guys, I was comparing my answer to something wolfram alpha gave me and it seemed different. I wish that I could accept one of the above commenters answer. – user1610950 Aug 01 '15 at 07:03
  • 1
    In the general case transposing a matrix *is not* inverting it (only for the special case of orthogonal matrices transposition is inversion). – datenwolf Aug 01 '15 at 09:39
  • 1
    see [matrix_inv in C++](http://stackoverflow.com/a/28084380/2521214) if you want to check the result just multiply original and inverted matrix and the result should be unit matrix ... (matrix_mul) – Spektre Aug 03 '15 at 05:52
  • @Spektre that's a great tip, thanks! – user1610950 Aug 03 '15 at 08:50

0 Answers0