当前位置： 文江博客话题详情

反转 4x4 矩阵

发布于 2024-07-27 19:37:13 字数 127 浏览 4 评论 0原文

我正在寻找有关如何反转 4x4 矩阵的示例代码实现。 我知道有高斯消除、LU 分解等，但我并没有详细研究它们，而是只是在寻找执行此操作的代码。

理想的语言是 C++，数据以列优先顺序以 16 个浮点数的数组形式提供。

原文

I am looking for a sample code implementation on how to invert a 4x4 matrix. I know there is Gaussian eleminiation, LU decomposition, etc., but instead of looking at them in detail I am really just looking for the code to do this.

Language ideally C++, data is available in array of 16 floats in column-major order.

如果你对这篇内容有疑问，欢迎到本站社区发帖提问 参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

扫码二维码加入Web技术交流群

发布评论

需要 登录 才能够评论， 你可以免费 注册 一个本站的账号。

评论（11）

平生欢 2024-08-03 19:37:13

此处：

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;
}

这是从 GLU 库的 MESA 实现中提取的。

here:

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;
}

This was lifted from MESA implementation of the GLU library.

回复 原文
傾城如夢未必闌珊 2024-08-03 19:37:13

如果有人寻找更定制的代码和“更易于阅读”，那么我得到了这个，

var A2323 = m.m22 * m.m33 - m.m23 * m.m32 ;
var A1323 = m.m21 * m.m33 - m.m23 * m.m31 ;
var A1223 = m.m21 * m.m32 - m.m22 * m.m31 ;
var A0323 = m.m20 * m.m33 - m.m23 * m.m30 ;
var A0223 = m.m20 * m.m32 - m.m22 * m.m30 ;
var A0123 = m.m20 * m.m31 - m.m21 * m.m30 ;
var A2313 = m.m12 * m.m33 - m.m13 * m.m32 ;
var A1313 = m.m11 * m.m33 - m.m13 * m.m31 ;
var A1213 = m.m11 * m.m32 - m.m12 * m.m31 ;
var A2312 = m.m12 * m.m23 - m.m13 * m.m22 ;
var A1312 = m.m11 * m.m23 - m.m13 * m.m21 ;
var A1212 = m.m11 * m.m22 - m.m12 * m.m21 ;
var A0313 = m.m10 * m.m33 - m.m13 * m.m30 ;
var A0213 = m.m10 * m.m32 - m.m12 * m.m30 ;
var A0312 = m.m10 * m.m23 - m.m13 * m.m20 ;
var A0212 = m.m10 * m.m22 - m.m12 * m.m20 ;
var A0113 = m.m10 * m.m31 - m.m11 * m.m30 ;
var A0112 = m.m10 * m.m21 - m.m11 * m.m20 ;

var det = m.m00 * ( m.m11 * A2323 - m.m12 * A1323 + m.m13 * A1223 ) 
    - m.m01 * ( m.m10 * A2323 - m.m12 * A0323 + m.m13 * A0223 ) 
    + m.m02 * ( m.m10 * A1323 - m.m11 * A0323 + m.m13 * A0123 ) 
    - m.m03 * ( m.m10 * A1223 - m.m11 * A0223 + m.m12 * A0123 ) ;
det = 1 / det;

return new Matrix4x4() {
   m00 = det *   ( m.m11 * A2323 - m.m12 * A1323 + m.m13 * A1223 ),
   m01 = det * - ( m.m01 * A2323 - m.m02 * A1323 + m.m03 * A1223 ),
   m02 = det *   ( m.m01 * A2313 - m.m02 * A1313 + m.m03 * A1213 ),
   m03 = det * - ( m.m01 * A2312 - m.m02 * A1312 + m.m03 * A1212 ),
   m10 = det * - ( m.m10 * A2323 - m.m12 * A0323 + m.m13 * A0223 ),
   m11 = det *   ( m.m00 * A2323 - m.m02 * A0323 + m.m03 * A0223 ),
   m12 = det * - ( m.m00 * A2313 - m.m02 * A0313 + m.m03 * A0213 ),
   m13 = det *   ( m.m00 * A2312 - m.m02 * A0312 + m.m03 * A0212 ),
   m20 = det *   ( m.m10 * A1323 - m.m11 * A0323 + m.m13 * A0123 ),
   m21 = det * - ( m.m00 * A1323 - m.m01 * A0323 + m.m03 * A0123 ),
   m22 = det *   ( m.m00 * A1313 - m.m01 * A0313 + m.m03 * A0113 ),
   m23 = det * - ( m.m00 * A1312 - m.m01 * A0312 + m.m03 * A0112 ),
   m30 = det * - ( m.m10 * A1223 - m.m11 * A0223 + m.m12 * A0123 ),
   m31 = det *   ( m.m00 * A1223 - m.m01 * A0223 + m.m02 * A0123 ),
   m32 = det * - ( m.m00 * A1213 - m.m01 * A0213 + m.m02 * A0113 ),
   m33 = det *   ( m.m00 * A1212 - m.m01 * A0212 + m.m02 * A0112 ),
};

我不写代码，但我的程序写了。 我编写了一个小程序来编写一个计算任意 N 矩阵的行列式和逆矩阵的程序。

我这样做是因为过去我需要一个反转 5x5 矩阵的代码，但地球上没有人这样做过，所以我做了一个。

此处查看该程序。

编辑：矩阵布局是逐行的（意味着 m01 位于第一行第二列）。 语言也是 C#，但应该很容易转换为 C。

If anyone looking for more costumized code and "easier to read", then I got this

var A2323 = m.m22 * m.m33 - m.m23 * m.m32 ;
var A1323 = m.m21 * m.m33 - m.m23 * m.m31 ;
var A1223 = m.m21 * m.m32 - m.m22 * m.m31 ;
var A0323 = m.m20 * m.m33 - m.m23 * m.m30 ;
var A0223 = m.m20 * m.m32 - m.m22 * m.m30 ;
var A0123 = m.m20 * m.m31 - m.m21 * m.m30 ;
var A2313 = m.m12 * m.m33 - m.m13 * m.m32 ;
var A1313 = m.m11 * m.m33 - m.m13 * m.m31 ;
var A1213 = m.m11 * m.m32 - m.m12 * m.m31 ;
var A2312 = m.m12 * m.m23 - m.m13 * m.m22 ;
var A1312 = m.m11 * m.m23 - m.m13 * m.m21 ;
var A1212 = m.m11 * m.m22 - m.m12 * m.m21 ;
var A0313 = m.m10 * m.m33 - m.m13 * m.m30 ;
var A0213 = m.m10 * m.m32 - m.m12 * m.m30 ;
var A0312 = m.m10 * m.m23 - m.m13 * m.m20 ;
var A0212 = m.m10 * m.m22 - m.m12 * m.m20 ;
var A0113 = m.m10 * m.m31 - m.m11 * m.m30 ;
var A0112 = m.m10 * m.m21 - m.m11 * m.m20 ;

var det = m.m00 * ( m.m11 * A2323 - m.m12 * A1323 + m.m13 * A1223 ) 
    - m.m01 * ( m.m10 * A2323 - m.m12 * A0323 + m.m13 * A0223 ) 
    + m.m02 * ( m.m10 * A1323 - m.m11 * A0323 + m.m13 * A0123 ) 
    - m.m03 * ( m.m10 * A1223 - m.m11 * A0223 + m.m12 * A0123 ) ;
det = 1 / det;

return new Matrix4x4() {
   m00 = det *   ( m.m11 * A2323 - m.m12 * A1323 + m.m13 * A1223 ),
   m01 = det * - ( m.m01 * A2323 - m.m02 * A1323 + m.m03 * A1223 ),
   m02 = det *   ( m.m01 * A2313 - m.m02 * A1313 + m.m03 * A1213 ),
   m03 = det * - ( m.m01 * A2312 - m.m02 * A1312 + m.m03 * A1212 ),
   m10 = det * - ( m.m10 * A2323 - m.m12 * A0323 + m.m13 * A0223 ),
   m11 = det *   ( m.m00 * A2323 - m.m02 * A0323 + m.m03 * A0223 ),
   m12 = det * - ( m.m00 * A2313 - m.m02 * A0313 + m.m03 * A0213 ),
   m13 = det *   ( m.m00 * A2312 - m.m02 * A0312 + m.m03 * A0212 ),
   m20 = det *   ( m.m10 * A1323 - m.m11 * A0323 + m.m13 * A0123 ),
   m21 = det * - ( m.m00 * A1323 - m.m01 * A0323 + m.m03 * A0123 ),
   m22 = det *   ( m.m00 * A1313 - m.m01 * A0313 + m.m03 * A0113 ),
   m23 = det * - ( m.m00 * A1312 - m.m01 * A0312 + m.m03 * A0112 ),
   m30 = det * - ( m.m10 * A1223 - m.m11 * A0223 + m.m12 * A0123 ),
   m31 = det *   ( m.m00 * A1223 - m.m01 * A0223 + m.m02 * A0123 ),
   m32 = det * - ( m.m00 * A1213 - m.m01 * A0213 + m.m02 * A0113 ),
   m33 = det *   ( m.m00 * A1212 - m.m01 * A0212 + m.m02 * A0112 ),
};

I don't write the code, but my program did. I made a small program to make a program that calculate the determinant and inverse of any N-matrix.

I do it because once in the past I need a code that inverses 5x5 matrix, but nobody in the earth have done this so I made one.

Take a look about the program here.

EDIT: The matrix layout is row-by-row (meaning m01 is in the first row and second column). Also the language is C#, but should be easy to convert into C.

回复 原文
多谢你的绝情让我学会死心 2024-08-03 19:37:13

我“汇总”了 MESA 实现（还编写了几个单元测试以确保它确实有效）。

在这里：

float invf(int i,int j,const float* m){

    int o = 2+(j-i);

    i += 4+o;
    j += 4-o;

    #define e(a,b) m[ ((j+b)%4)*4 + ((i+a)%4) ]

    float inv =
     + e(+1,-1)*e(+0,+0)*e(-1,+1)
     + e(+1,+1)*e(+0,-1)*e(-1,+0)
     + e(-1,-1)*e(+1,+0)*e(+0,+1)
     - e(-1,-1)*e(+0,+0)*e(+1,+1)
     - e(-1,+1)*e(+0,-1)*e(+1,+0)
     - e(+1,-1)*e(-1,+0)*e(+0,+1);

    return (o%2)?inv : -inv;

    #undef e

}

bool inverseMatrix4x4(const float *m, float *out)
{

    float inv[16];

    for(int i=0;i<4;i++)
        for(int j=0;j<4;j++)
            inv[j*4+i] = invf(i,j,m);

    double D = 0;

    for(int k=0;k<4;k++) D += m[k] * inv[k*4];

    if (D == 0) return false;

    D = 1.0 / D;

    for (int i = 0; i < 16; i++)
        out[i] = inv[i] * D;

    return true;

}

我写了一些关于此的内容，并显示了正/负因素的模式 在我的博客上。

正如@LiraNuna 所建议的，在许多平台上都可以使用此类例程的硬件加速版本，因此我很高兴拥有一个可读且简洁的“备份版本”。

注意：这可能比 MESA 实现慢 3.5 倍或更差。 您可以改变因素的模式以删除一些添加内容等...但它会降低可读性并且仍然不会很快。

I 'rolled up' the MESA implementation (also wrote a couple of unit tests to ensure it actually works).

Here:

float invf(int i,int j,const float* m){

    int o = 2+(j-i);

    i += 4+o;
    j += 4-o;

    #define e(a,b) m[ ((j+b)%4)*4 + ((i+a)%4) ]

    float inv =
     + e(+1,-1)*e(+0,+0)*e(-1,+1)
     + e(+1,+1)*e(+0,-1)*e(-1,+0)
     + e(-1,-1)*e(+1,+0)*e(+0,+1)
     - e(-1,-1)*e(+0,+0)*e(+1,+1)
     - e(-1,+1)*e(+0,-1)*e(+1,+0)
     - e(+1,-1)*e(-1,+0)*e(+0,+1);

    return (o%2)?inv : -inv;

    #undef e

}

bool inverseMatrix4x4(const float *m, float *out)
{

    float inv[16];

    for(int i=0;i<4;i++)
        for(int j=0;j<4;j++)
            inv[j*4+i] = invf(i,j,m);

    double D = 0;

    for(int k=0;k<4;k++) D += m[k] * inv[k*4];

    if (D == 0) return false;

    D = 1.0 / D;

    for (int i = 0; i < 16; i++)
        out[i] = inv[i] * D;

    return true;

}

I wrote a little about this and display the pattern of positive/negative factors on my blog.

As suggested by @LiraNuna, on many platforms hardware accelerated versions of such routines are available so I'm happy to have a 'backup version' that's readable and concise.

Note: this may run 3.5 times slower or worse than the MESA implementation. You can shift the pattern of factors to remove some additions etc... but it would lose in readability and still won't be very fast.

回复 原文
小忆控 2024-08-03 19:37:13

如果您需要具有大量函数的 C++ 矩阵库，请查看 Eigen 库 - http://eigen.tuxfamily.org

If you need a C++ matrix library with a lot of functions, have a look at Eigen library - http://eigen.tuxfamily.org

回复 原文
痕至 2024-08-03 19:37:13

这是 @willnode 的答案的 C++ 版本

template<typename Matrix>
static inline void InvertMatrix4(const Matrix& m, Matrix& im, double& det)
{
    double A2323 = m(2, 2) * m(3, 3) - m(2, 3) * m(3, 2);
    double A1323 = m(2, 1) * m(3, 3) - m(2, 3) * m(3, 1);
    double A1223 = m(2, 1) * m(3, 2) - m(2, 2) * m(3, 1);
    double A0323 = m(2, 0) * m(3, 3) - m(2, 3) * m(3, 0);
    double A0223 = m(2, 0) * m(3, 2) - m(2, 2) * m(3, 0);
    double A0123 = m(2, 0) * m(3, 1) - m(2, 1) * m(3, 0);
    double A2313 = m(1, 2) * m(3, 3) - m(1, 3) * m(3, 2);
    double A1313 = m(1, 1) * m(3, 3) - m(1, 3) * m(3, 1);
    double A1213 = m(1, 1) * m(3, 2) - m(1, 2) * m(3, 1);
    double A2312 = m(1, 2) * m(2, 3) - m(1, 3) * m(2, 2);
    double A1312 = m(1, 1) * m(2, 3) - m(1, 3) * m(2, 1);
    double A1212 = m(1, 1) * m(2, 2) - m(1, 2) * m(2, 1);
    double A0313 = m(1, 0) * m(3, 3) - m(1, 3) * m(3, 0);
    double A0213 = m(1, 0) * m(3, 2) - m(1, 2) * m(3, 0);
    double A0312 = m(1, 0) * m(2, 3) - m(1, 3) * m(2, 0);
    double A0212 = m(1, 0) * m(2, 2) - m(1, 2) * m(2, 0);
    double A0113 = m(1, 0) * m(3, 1) - m(1, 1) * m(3, 0);
    double A0112 = m(1, 0) * m(2, 1) - m(1, 1) * m(2, 0);

    det = m(0, 0) * ( m(1, 1) * A2323 - m(1, 2) * A1323 + m(1, 3) * A1223 )
        - m(0, 1) * ( m(1, 0) * A2323 - m(1, 2) * A0323 + m(1, 3) * A0223 )
        + m(0, 2) * ( m(1, 0) * A1323 - m(1, 1) * A0323 + m(1, 3) * A0123 )
        - m(0, 3) * ( m(1, 0) * A1223 - m(1, 1) * A0223 + m(1, 2) * A0123 );
    det = 1 / det;

    im(0, 0) = det *   ( m(1, 1) * A2323 - m(1, 2) * A1323 + m(1, 3) * A1223 );
    im(0, 1) = det * - ( m(0, 1) * A2323 - m(0, 2) * A1323 + m(0, 3) * A1223 );
    im(0, 2) = det *   ( m(0, 1) * A2313 - m(0, 2) * A1313 + m(0, 3) * A1213 );
    im(0, 3) = det * - ( m(0, 1) * A2312 - m(0, 2) * A1312 + m(0, 3) * A1212 );
    im(1, 0) = det * - ( m(1, 0) * A2323 - m(1, 2) * A0323 + m(1, 3) * A0223 );
    im(1, 1) = det *   ( m(0, 0) * A2323 - m(0, 2) * A0323 + m(0, 3) * A0223 );
    im(1, 2) = det * - ( m(0, 0) * A2313 - m(0, 2) * A0313 + m(0, 3) * A0213 );
    im(1, 3) = det *   ( m(0, 0) * A2312 - m(0, 2) * A0312 + m(0, 3) * A0212 );
    im(2, 0) = det *   ( m(1, 0) * A1323 - m(1, 1) * A0323 + m(1, 3) * A0123 );
    im(2, 1) = det * - ( m(0, 0) * A1323 - m(0, 1) * A0323 + m(0, 3) * A0123 );
    im(2, 2) = det *   ( m(0, 0) * A1313 - m(0, 1) * A0313 + m(0, 3) * A0113 );
    im(2, 3) = det * - ( m(0, 0) * A1312 - m(0, 1) * A0312 + m(0, 3) * A0112 );
    im(3, 0) = det * - ( m(1, 0) * A1223 - m(1, 1) * A0223 + m(1, 2) * A0123 );
    im(3, 1) = det *   ( m(0, 0) * A1223 - m(0, 1) * A0223 + m(0, 2) * A0123 );
    im(3, 2) = det * - ( m(0, 0) * A1213 - m(0, 1) * A0213 + m(0, 2) * A0113 );
    im(3, 3) = det *   ( m(0, 0) * A1212 - m(0, 1) * A0212 + m(0, 2) * A0112 );
}

This is the C++ version for @willnode's answer

template<typename Matrix>
static inline void InvertMatrix4(const Matrix& m, Matrix& im, double& det)
{
    double A2323 = m(2, 2) * m(3, 3) - m(2, 3) * m(3, 2);
    double A1323 = m(2, 1) * m(3, 3) - m(2, 3) * m(3, 1);
    double A1223 = m(2, 1) * m(3, 2) - m(2, 2) * m(3, 1);
    double A0323 = m(2, 0) * m(3, 3) - m(2, 3) * m(3, 0);
    double A0223 = m(2, 0) * m(3, 2) - m(2, 2) * m(3, 0);
    double A0123 = m(2, 0) * m(3, 1) - m(2, 1) * m(3, 0);
    double A2313 = m(1, 2) * m(3, 3) - m(1, 3) * m(3, 2);
    double A1313 = m(1, 1) * m(3, 3) - m(1, 3) * m(3, 1);
    double A1213 = m(1, 1) * m(3, 2) - m(1, 2) * m(3, 1);
    double A2312 = m(1, 2) * m(2, 3) - m(1, 3) * m(2, 2);
    double A1312 = m(1, 1) * m(2, 3) - m(1, 3) * m(2, 1);
    double A1212 = m(1, 1) * m(2, 2) - m(1, 2) * m(2, 1);
    double A0313 = m(1, 0) * m(3, 3) - m(1, 3) * m(3, 0);
    double A0213 = m(1, 0) * m(3, 2) - m(1, 2) * m(3, 0);
    double A0312 = m(1, 0) * m(2, 3) - m(1, 3) * m(2, 0);
    double A0212 = m(1, 0) * m(2, 2) - m(1, 2) * m(2, 0);
    double A0113 = m(1, 0) * m(3, 1) - m(1, 1) * m(3, 0);
    double A0112 = m(1, 0) * m(2, 1) - m(1, 1) * m(2, 0);

    det = m(0, 0) * ( m(1, 1) * A2323 - m(1, 2) * A1323 + m(1, 3) * A1223 )
        - m(0, 1) * ( m(1, 0) * A2323 - m(1, 2) * A0323 + m(1, 3) * A0223 )
        + m(0, 2) * ( m(1, 0) * A1323 - m(1, 1) * A0323 + m(1, 3) * A0123 )
        - m(0, 3) * ( m(1, 0) * A1223 - m(1, 1) * A0223 + m(1, 2) * A0123 );
    det = 1 / det;

    im(0, 0) = det *   ( m(1, 1) * A2323 - m(1, 2) * A1323 + m(1, 3) * A1223 );
    im(0, 1) = det * - ( m(0, 1) * A2323 - m(0, 2) * A1323 + m(0, 3) * A1223 );
    im(0, 2) = det *   ( m(0, 1) * A2313 - m(0, 2) * A1313 + m(0, 3) * A1213 );
    im(0, 3) = det * - ( m(0, 1) * A2312 - m(0, 2) * A1312 + m(0, 3) * A1212 );
    im(1, 0) = det * - ( m(1, 0) * A2323 - m(1, 2) * A0323 + m(1, 3) * A0223 );
    im(1, 1) = det *   ( m(0, 0) * A2323 - m(0, 2) * A0323 + m(0, 3) * A0223 );
    im(1, 2) = det * - ( m(0, 0) * A2313 - m(0, 2) * A0313 + m(0, 3) * A0213 );
    im(1, 3) = det *   ( m(0, 0) * A2312 - m(0, 2) * A0312 + m(0, 3) * A0212 );
    im(2, 0) = det *   ( m(1, 0) * A1323 - m(1, 1) * A0323 + m(1, 3) * A0123 );
    im(2, 1) = det * - ( m(0, 0) * A1323 - m(0, 1) * A0323 + m(0, 3) * A0123 );
    im(2, 2) = det *   ( m(0, 0) * A1313 - m(0, 1) * A0313 + m(0, 3) * A0113 );
    im(2, 3) = det * - ( m(0, 0) * A1312 - m(0, 1) * A0312 + m(0, 3) * A0112 );
    im(3, 0) = det * - ( m(1, 0) * A1223 - m(1, 1) * A0223 + m(1, 2) * A0123 );
    im(3, 1) = det *   ( m(0, 0) * A1223 - m(0, 1) * A0223 + m(0, 2) * A0123 );
    im(3, 2) = det * - ( m(0, 0) * A1213 - m(0, 1) * A0213 + m(0, 2) * A0113 );
    im(3, 3) = det *   ( m(0, 0) * A1212 - m(0, 1) * A0212 + m(0, 2) * A0112 );
}
回复 原文
伤痕我心 2024-08-03 19:37:13

您可以使用 GNU 科学库 或在其中查找代码。

编辑：您似乎想要 线性代数 部分。

You can use the GNU Scientific Library or look the code up in it.

Edit: You seem to want the Linear Algebra section.

回复 原文
黒涩兲箜 2024-08-03 19:37:13

受 @shoosh 的启发，我检查了 MESA 实现，发现矩阵求逆在最近的 mesa 版本中看起来非常不同。 我想这些都是很好的改进。 这是来自 Mesa-17.3.9 的矩阵求逆代码：

/* Returns true for success, false for failure (singular matrix) */
bool DirectVolumeRenderer::_mesa_invert_matrix_general( GLfloat out[16], const GLfloat in[16] )
{
    /**
     * References an element of 4x4 matrix.
     * Calculate the linear storage index of the element and references it. 
     */
    #define MAT(m,r,c) (m)[(c)*4+(r)]
    /**
     * Swaps the values of two floating point variables.
     */
    #define SWAP_ROWS(a, b) { GLfloat *_tmp = a; (a)=(b); (b)=_tmp; }

    const GLfloat *m = in;
    GLfloat wtmp[4][8];
    GLfloat m0, m1, m2, m3, s;
    GLfloat *r0, *r1, *r2, *r3;

    r0 = wtmp[0], r1 = wtmp[1], r2 = wtmp[2], r3 = wtmp[3];

    r0[0] = MAT(m,0,0), r0[1] = MAT(m,0,1),
    r0[2] = MAT(m,0,2), r0[3] = MAT(m,0,3),
    r0[4] = 1.0, r0[5] = r0[6] = r0[7] = 0.0,

    r1[0] = MAT(m,1,0), r1[1] = MAT(m,1,1),
    r1[2] = MAT(m,1,2), r1[3] = MAT(m,1,3),
    r1[5] = 1.0, r1[4] = r1[6] = r1[7] = 0.0,

    r2[0] = MAT(m,2,0), r2[1] = MAT(m,2,1),
    r2[2] = MAT(m,2,2), r2[3] = MAT(m,2,3),
    r2[6] = 1.0, r2[4] = r2[5] = r2[7] = 0.0,

    r3[0] = MAT(m,3,0), r3[1] = MAT(m,3,1),
    r3[2] = MAT(m,3,2), r3[3] = MAT(m,3,3),
    r3[7] = 1.0, r3[4] = r3[5] = r3[6] = 0.0;

    /* choose pivot - or die */
    if (fabsf(r3[0])>fabsf(r2[0])) SWAP_ROWS(r3, r2);
    if (fabsf(r2[0])>fabsf(r1[0])) SWAP_ROWS(r2, r1);
    if (fabsf(r1[0])>fabsf(r0[0])) SWAP_ROWS(r1, r0);
    if (0.0F == r0[0])
        return false;

    /* eliminate first variable     */
    m1 = r1[0]/r0[0]; m2 = r2[0]/r0[0]; m3 = r3[0]/r0[0];
    s = r0[1]; r1[1] -= m1 * s; r2[1] -= m2 * s; r3[1] -= m3 * s;
    s = r0[2]; r1[2] -= m1 * s; r2[2] -= m2 * s; r3[2] -= m3 * s;
    s = r0[3]; r1[3] -= m1 * s; r2[3] -= m2 * s; r3[3] -= m3 * s;
    s = r0[4];
    if (s != 0.0F) { r1[4] -= m1 * s; r2[4] -= m2 * s; r3[4] -= m3 * s; }
    s = r0[5];
    if (s != 0.0F) { r1[5] -= m1 * s; r2[5] -= m2 * s; r3[5] -= m3 * s; }
    s = r0[6];
    if (s != 0.0F) { r1[6] -= m1 * s; r2[6] -= m2 * s; r3[6] -= m3 * s; }
    s = r0[7];
    if (s != 0.0F) { r1[7] -= m1 * s; r2[7] -= m2 * s; r3[7] -= m3 * s; }

    /* choose pivot - or die */
    if (fabsf(r3[1])>fabsf(r2[1])) SWAP_ROWS(r3, r2);
    if (fabsf(r2[1])>fabsf(r1[1])) SWAP_ROWS(r2, r1);
    if (0.0F == r1[1])
        return false;

    /* eliminate second variable */
    m2 = r2[1]/r1[1]; m3 = r3[1]/r1[1];
    r2[2] -= m2 * r1[2]; r3[2] -= m3 * r1[2];
    r2[3] -= m2 * r1[3]; r3[3] -= m3 * r1[3];
    s = r1[4]; if (0.0F != s) { r2[4] -= m2 * s; r3[4] -= m3 * s; }
    s = r1[5]; if (0.0F != s) { r2[5] -= m2 * s; r3[5] -= m3 * s; }
    s = r1[6]; if (0.0F != s) { r2[6] -= m2 * s; r3[6] -= m3 * s; }
    s = r1[7]; if (0.0F != s) { r2[7] -= m2 * s; r3[7] -= m3 * s; }

    /* choose pivot - or die */
    if (fabsf(r3[2])>fabsf(r2[2])) SWAP_ROWS(r3, r2);
    if (0.0F == r2[2])
        return false;

    /* eliminate third variable */
    m3 = r3[2]/r2[2];
    r3[3] -= m3 * r2[3], r3[4] -= m3 * r2[4],
    r3[5] -= m3 * r2[5], r3[6] -= m3 * r2[6],
    r3[7] -= m3 * r2[7];

    /* last check */
    if (0.0F == r3[3])
        return false;

    s = 1.0F/r3[3];             /* now back substitute row 3 */
    r3[4] *= s; r3[5] *= s; r3[6] *= s; r3[7] *= s;

    m2 = r2[3];                 /* now back substitute row 2 */
    s  = 1.0F/r2[2];
    r2[4] = s * (r2[4] - r3[4] * m2), r2[5] = s * (r2[5] - r3[5] * m2),
    r2[6] = s * (r2[6] - r3[6] * m2), r2[7] = s * (r2[7] - r3[7] * m2);
    m1 = r1[3];
    r1[4] -= r3[4] * m1, r1[5] -= r3[5] * m1,
    r1[6] -= r3[6] * m1, r1[7] -= r3[7] * m1;
    m0 = r0[3];
    r0[4] -= r3[4] * m0, r0[5] -= r3[5] * m0,
    r0[6] -= r3[6] * m0, r0[7] -= r3[7] * m0;

    m1 = r1[2];                 /* now back substitute row 1 */
    s  = 1.0F/r1[1];
    r1[4] = s * (r1[4] - r2[4] * m1), r1[5] = s * (r1[5] - r2[5] * m1),
    r1[6] = s * (r1[6] - r2[6] * m1), r1[7] = s * (r1[7] - r2[7] * m1);
    m0 = r0[2];
    r0[4] -= r2[4] * m0, r0[5] -= r2[5] * m0,
    r0[6] -= r2[6] * m0, r0[7] -= r2[7] * m0;

    m0 = r0[1];                 /* now back substitute row 0 */
    s  = 1.0F/r0[0];
    r0[4] = s * (r0[4] - r1[4] * m0), r0[5] = s * (r0[5] - r1[5] * m0),
    r0[6] = s * (r0[6] - r1[6] * m0), r0[7] = s * (r0[7] - r1[7] * m0);

    MAT(out,0,0) = r0[4]; MAT(out,0,1) = r0[5],
    MAT(out,0,2) = r0[6]; MAT(out,0,3) = r0[7],
    MAT(out,1,0) = r1[4]; MAT(out,1,1) = r1[5],
    MAT(out,1,2) = r1[6]; MAT(out,1,3) = r1[7],
    MAT(out,2,0) = r2[4]; MAT(out,2,1) = r2[5],
    MAT(out,2,2) = r2[6]; MAT(out,2,3) = r2[7],
    MAT(out,3,0) = r3[4]; MAT(out,3,1) = r3[5],
    MAT(out,3,2) = r3[6]; MAT(out,3,3) = r3[7];

    #undef SWAP_ROWS
    #undef MAT

    return true;
}

注意：您可以在 mesa 代码库中找到这段代码：mesa-17.3.9/src/mesa/math/m_matrix.c。

Inspired by @shoosh to check out MESA implementations, I found that matrix inversion looks quite different in more recent mesa releases. I suppose those are good improvements. Here's the matrix inversion code from Mesa-17.3.9:

/* Returns true for success, false for failure (singular matrix) */
bool DirectVolumeRenderer::_mesa_invert_matrix_general( GLfloat out[16], const GLfloat in[16] )
{
    /**
     * References an element of 4x4 matrix.
     * Calculate the linear storage index of the element and references it. 
     */
    #define MAT(m,r,c) (m)[(c)*4+(r)]
    /**
     * Swaps the values of two floating point variables.
     */
    #define SWAP_ROWS(a, b) { GLfloat *_tmp = a; (a)=(b); (b)=_tmp; }

    const GLfloat *m = in;
    GLfloat wtmp[4][8];
    GLfloat m0, m1, m2, m3, s;
    GLfloat *r0, *r1, *r2, *r3;

    r0 = wtmp[0], r1 = wtmp[1], r2 = wtmp[2], r3 = wtmp[3];

    r0[0] = MAT(m,0,0), r0[1] = MAT(m,0,1),
    r0[2] = MAT(m,0,2), r0[3] = MAT(m,0,3),
    r0[4] = 1.0, r0[5] = r0[6] = r0[7] = 0.0,

    r1[0] = MAT(m,1,0), r1[1] = MAT(m,1,1),
    r1[2] = MAT(m,1,2), r1[3] = MAT(m,1,3),
    r1[5] = 1.0, r1[4] = r1[6] = r1[7] = 0.0,

    r2[0] = MAT(m,2,0), r2[1] = MAT(m,2,1),
    r2[2] = MAT(m,2,2), r2[3] = MAT(m,2,3),
    r2[6] = 1.0, r2[4] = r2[5] = r2[7] = 0.0,

    r3[0] = MAT(m,3,0), r3[1] = MAT(m,3,1),
    r3[2] = MAT(m,3,2), r3[3] = MAT(m,3,3),
    r3[7] = 1.0, r3[4] = r3[5] = r3[6] = 0.0;

    /* choose pivot - or die */
    if (fabsf(r3[0])>fabsf(r2[0])) SWAP_ROWS(r3, r2);
    if (fabsf(r2[0])>fabsf(r1[0])) SWAP_ROWS(r2, r1);
    if (fabsf(r1[0])>fabsf(r0[0])) SWAP_ROWS(r1, r0);
    if (0.0F == r0[0])
        return false;

    /* eliminate first variable     */
    m1 = r1[0]/r0[0]; m2 = r2[0]/r0[0]; m3 = r3[0]/r0[0];
    s = r0[1]; r1[1] -= m1 * s; r2[1] -= m2 * s; r3[1] -= m3 * s;
    s = r0[2]; r1[2] -= m1 * s; r2[2] -= m2 * s; r3[2] -= m3 * s;
    s = r0[3]; r1[3] -= m1 * s; r2[3] -= m2 * s; r3[3] -= m3 * s;
    s = r0[4];
    if (s != 0.0F) { r1[4] -= m1 * s; r2[4] -= m2 * s; r3[4] -= m3 * s; }
    s = r0[5];
    if (s != 0.0F) { r1[5] -= m1 * s; r2[5] -= m2 * s; r3[5] -= m3 * s; }
    s = r0[6];
    if (s != 0.0F) { r1[6] -= m1 * s; r2[6] -= m2 * s; r3[6] -= m3 * s; }
    s = r0[7];
    if (s != 0.0F) { r1[7] -= m1 * s; r2[7] -= m2 * s; r3[7] -= m3 * s; }

    /* choose pivot - or die */
    if (fabsf(r3[1])>fabsf(r2[1])) SWAP_ROWS(r3, r2);
    if (fabsf(r2[1])>fabsf(r1[1])) SWAP_ROWS(r2, r1);
    if (0.0F == r1[1])
        return false;

    /* eliminate second variable */
    m2 = r2[1]/r1[1]; m3 = r3[1]/r1[1];
    r2[2] -= m2 * r1[2]; r3[2] -= m3 * r1[2];
    r2[3] -= m2 * r1[3]; r3[3] -= m3 * r1[3];
    s = r1[4]; if (0.0F != s) { r2[4] -= m2 * s; r3[4] -= m3 * s; }
    s = r1[5]; if (0.0F != s) { r2[5] -= m2 * s; r3[5] -= m3 * s; }
    s = r1[6]; if (0.0F != s) { r2[6] -= m2 * s; r3[6] -= m3 * s; }
    s = r1[7]; if (0.0F != s) { r2[7] -= m2 * s; r3[7] -= m3 * s; }

    /* choose pivot - or die */
    if (fabsf(r3[2])>fabsf(r2[2])) SWAP_ROWS(r3, r2);
    if (0.0F == r2[2])
        return false;

    /* eliminate third variable */
    m3 = r3[2]/r2[2];
    r3[3] -= m3 * r2[3], r3[4] -= m3 * r2[4],
    r3[5] -= m3 * r2[5], r3[6] -= m3 * r2[6],
    r3[7] -= m3 * r2[7];

    /* last check */
    if (0.0F == r3[3])
        return false;

    s = 1.0F/r3[3];             /* now back substitute row 3 */
    r3[4] *= s; r3[5] *= s; r3[6] *= s; r3[7] *= s;

    m2 = r2[3];                 /* now back substitute row 2 */
    s  = 1.0F/r2[2];
    r2[4] = s * (r2[4] - r3[4] * m2), r2[5] = s * (r2[5] - r3[5] * m2),
    r2[6] = s * (r2[6] - r3[6] * m2), r2[7] = s * (r2[7] - r3[7] * m2);
    m1 = r1[3];
    r1[4] -= r3[4] * m1, r1[5] -= r3[5] * m1,
    r1[6] -= r3[6] * m1, r1[7] -= r3[7] * m1;
    m0 = r0[3];
    r0[4] -= r3[4] * m0, r0[5] -= r3[5] * m0,
    r0[6] -= r3[6] * m0, r0[7] -= r3[7] * m0;

    m1 = r1[2];                 /* now back substitute row 1 */
    s  = 1.0F/r1[1];
    r1[4] = s * (r1[4] - r2[4] * m1), r1[5] = s * (r1[5] - r2[5] * m1),
    r1[6] = s * (r1[6] - r2[6] * m1), r1[7] = s * (r1[7] - r2[7] * m1);
    m0 = r0[2];
    r0[4] -= r2[4] * m0, r0[5] -= r2[5] * m0,
    r0[6] -= r2[6] * m0, r0[7] -= r2[7] * m0;

    m0 = r0[1];                 /* now back substitute row 0 */
    s  = 1.0F/r0[0];
    r0[4] = s * (r0[4] - r1[4] * m0), r0[5] = s * (r0[5] - r1[5] * m0),
    r0[6] = s * (r0[6] - r1[6] * m0), r0[7] = s * (r0[7] - r1[7] * m0);

    MAT(out,0,0) = r0[4]; MAT(out,0,1) = r0[5],
    MAT(out,0,2) = r0[6]; MAT(out,0,3) = r0[7],
    MAT(out,1,0) = r1[4]; MAT(out,1,1) = r1[5],
    MAT(out,1,2) = r1[6]; MAT(out,1,3) = r1[7],
    MAT(out,2,0) = r2[4]; MAT(out,2,1) = r2[5],
    MAT(out,2,2) = r2[6]; MAT(out,2,3) = r2[7],
    MAT(out,3,0) = r3[4]; MAT(out,3,1) = r3[5],
    MAT(out,3,2) = r3[6]; MAT(out,3,3) = r3[7];

    #undef SWAP_ROWS
    #undef MAT

    return true;
}

Note: you can find this piece of code in the mesa code base: mesa-17.3.9/src/mesa/math/m_matrix.c.

回复 原文
夏天碎花小短裙 2024-08-03 19:37:13

这是一个小型（只有一个标头）C++ 向量数学 库（齿轮面向 3D 编程）。 如果你使用它，请记住，与 OpenGL 期望的相比，它在内存中的矩阵布局是相反的，我很高兴弄清楚它......

Here is a small (just one header) C++ vector math library (geared towards 3D programming). If you use it, keep in mind that layout of its matrices in memory is inverted comparing to what OpenGL expects, I had fun time figuring it out...

回复 原文
风轻花落早 2024-08-03 19:37:13

您可以根据此博客使其更快。

#define SUBP(i,j) input[i][j]
#define SUBQ(i,j) input[i][2+j]
#define SUBR(i,j) input[2+i][j]
#define SUBS(i,j) input[2+i][2+j]

#define OUTP(i,j) output[i][j]
#define OUTQ(i,j) output[i][2+j]
#define OUTR(i,j) output[2+i][j]
#define OUTS(i,j) output[2+i][2+j]

#define INVP(i,j) invP[i][j]
#define INVPQ(i,j) invPQ[i][j]
#define RINVP(i,j) RinvP[i][j]
#define INVPQ(i,j) invPQ[i][j]
#define RINVPQ(i,j) RinvPQ[i][j]
#define INVPQR(i,j) invPQR[i][j]
#define INVS(i,j) invS[i][j]

#define MULTI(MAT1, MAT2, MAT3) \
    MAT3(0,0)=MAT1(0,0)*MAT2(0,0) + MAT1(0,1)*MAT2(1,0); \
MAT3(0,1)=MAT1(0,0)*MAT2(0,1) + MAT1(0,1)*MAT2(1,1); \
MAT3(1,0)=MAT1(1,0)*MAT2(0,0) + MAT1(1,1)*MAT2(1,0); \
MAT3(1,1)=MAT1(1,0)*MAT2(0,1) + MAT1(1,1)*MAT2(1,1);

#define INV(MAT1, MAT2) \
    _det = 1.0 / (MAT1(0,0) * MAT1(1,1) - MAT1(0,1) * MAT1(1,0)); \
MAT2(0,0) = MAT1(1,1) * _det; \
MAT2(1,1) = MAT1(0,0) * _det; \
MAT2(0,1) = -MAT1(0,1) * _det; \
MAT2(1,0) = -MAT1(1,0) * _det; \

#define SUBTRACT(MAT1, MAT2, MAT3) \
    MAT3(0,0)=MAT1(0,0) - MAT2(0,0); \
MAT3(0,1)=MAT1(0,1) - MAT2(0,1); \
MAT3(1,0)=MAT1(1,0) - MAT2(1,0); \
MAT3(1,1)=MAT1(1,1) - MAT2(1,1);

#define NEGATIVE(MAT) \
    MAT(0,0)=-MAT(0,0); \
MAT(0,1)=-MAT(0,1); \
MAT(1,0)=-MAT(1,0); \
MAT(1,1)=-MAT(1,1);


void getInvertMatrix(complex<double> input[4][4], complex<double> output[4][4]) {
    complex<double> _det;
    complex<double> invP[2][2];
    complex<double> invPQ[2][2];
    complex<double> RinvP[2][2];
    complex<double> RinvPQ[2][2];
    complex<double> invPQR[2][2];
    complex<double> invS[2][2];


    INV(SUBP, INVP);
    MULTI(SUBR, INVP, RINVP);
    MULTI(INVP, SUBQ, INVPQ);
    MULTI(RINVP, SUBQ, RINVPQ);
    SUBTRACT(SUBS, RINVPQ, INVS);
    INV(INVS, OUTS);
    NEGATIVE(OUTS);
    MULTI(OUTS, RINVP, OUTR);
    MULTI(INVPQ, OUTS, OUTQ);
    MULTI(INVPQ, OUTR, INVPQR);
    SUBTRACT(INVP, INVPQR, OUTP);
}

这不是一个完整的实现，因为 P 可能不可逆，但您可以将此代码与 MESA 实现结合起来以获得更好的性能。

You can make it faster according to this blog.

#define SUBP(i,j) input[i][j]
#define SUBQ(i,j) input[i][2+j]
#define SUBR(i,j) input[2+i][j]
#define SUBS(i,j) input[2+i][2+j]

#define OUTP(i,j) output[i][j]
#define OUTQ(i,j) output[i][2+j]
#define OUTR(i,j) output[2+i][j]
#define OUTS(i,j) output[2+i][2+j]

#define INVP(i,j) invP[i][j]
#define INVPQ(i,j) invPQ[i][j]
#define RINVP(i,j) RinvP[i][j]
#define INVPQ(i,j) invPQ[i][j]
#define RINVPQ(i,j) RinvPQ[i][j]
#define INVPQR(i,j) invPQR[i][j]
#define INVS(i,j) invS[i][j]

#define MULTI(MAT1, MAT2, MAT3) \
    MAT3(0,0)=MAT1(0,0)*MAT2(0,0) + MAT1(0,1)*MAT2(1,0); \
MAT3(0,1)=MAT1(0,0)*MAT2(0,1) + MAT1(0,1)*MAT2(1,1); \
MAT3(1,0)=MAT1(1,0)*MAT2(0,0) + MAT1(1,1)*MAT2(1,0); \
MAT3(1,1)=MAT1(1,0)*MAT2(0,1) + MAT1(1,1)*MAT2(1,1);

#define INV(MAT1, MAT2) \
    _det = 1.0 / (MAT1(0,0) * MAT1(1,1) - MAT1(0,1) * MAT1(1,0)); \
MAT2(0,0) = MAT1(1,1) * _det; \
MAT2(1,1) = MAT1(0,0) * _det; \
MAT2(0,1) = -MAT1(0,1) * _det; \
MAT2(1,0) = -MAT1(1,0) * _det; \

#define SUBTRACT(MAT1, MAT2, MAT3) \
    MAT3(0,0)=MAT1(0,0) - MAT2(0,0); \
MAT3(0,1)=MAT1(0,1) - MAT2(0,1); \
MAT3(1,0)=MAT1(1,0) - MAT2(1,0); \
MAT3(1,1)=MAT1(1,1) - MAT2(1,1);

#define NEGATIVE(MAT) \
    MAT(0,0)=-MAT(0,0); \
MAT(0,1)=-MAT(0,1); \
MAT(1,0)=-MAT(1,0); \
MAT(1,1)=-MAT(1,1);


void getInvertMatrix(complex<double> input[4][4], complex<double> output[4][4]) {
    complex<double> _det;
    complex<double> invP[2][2];
    complex<double> invPQ[2][2];
    complex<double> RinvP[2][2];
    complex<double> RinvPQ[2][2];
    complex<double> invPQR[2][2];
    complex<double> invS[2][2];


    INV(SUBP, INVP);
    MULTI(SUBR, INVP, RINVP);
    MULTI(INVP, SUBQ, INVPQ);
    MULTI(RINVP, SUBQ, RINVPQ);
    SUBTRACT(SUBS, RINVPQ, INVS);
    INV(INVS, OUTS);
    NEGATIVE(OUTS);
    MULTI(OUTS, RINVP, OUTR);
    MULTI(INVPQ, OUTS, OUTQ);
    MULTI(INVPQ, OUTR, INVPQR);
    SUBTRACT(INVP, INVPQR, OUTP);
}

This is not a complete implementation because P may not be invertible, but you can combine this code with MESA implementation to get a better performance.

回复 原文
雅心素梦 2024-08-03 19:37:13

添加一个可能对其他人有用的二维案例：

inline bool invert4by4matrix(mat m, mat inv_m)
 {

        inv_m->v[0][0] = m->v[1][1] * m->v[2][2] * m->v[3][3] -
                         m->v[1][1] * m->v[2][3] * m->v[3][2] -
                         m->v[2][1] * m->v[1][2] * m->v[3][3] +
                         m->v[2][1] * m->v[1][3] * m->v[3][2] +
                         m->v[3][1] * m->v[1][2] * m->v[2][3] -
                         m->v[3][1] * m->v[1][3] * m->v[2][2];

        inv_m->v[1][0] = -m->v[1][0] * m->v[2][2] * m->v[3][3] +
                         m->v[1][0] * m->v[2][3] * m->v[3][2] +
                         m->v[2][0] * m->v[1][2] * m->v[3][3] -
                         m->v[2][0] * m->v[1][3] * m->v[3][2] -
                         m->v[3][0] * m->v[1][2] * m->v[2][3] +
                         m->v[3][0] * m->v[1][3] * m->v[2][2];

        inv_m->v[2][0] = m->v[1][0] * m->v[2][1] * m->v[3][3] -
                         m->v[1][0] * m->v[2][3] * m->v[3][1] -
                         m->v[2][0] * m->v[1][1] * m->v[3][3] +
                         m->v[2][0] * m->v[1][3] * m->v[3][1] +
                         m->v[3][0] * m->v[1][1] * m->v[2][3] -
                         m->v[3][0] * m->v[1][3] * m->v[2][1];

        inv_m->v[3][0] = -m->v[1][0] * m->v[2][1] * m->v[3][2] +
                         m->v[1][0] * m->v[2][2] * m->v[3][1] +
                         m->v[2][0] * m->v[1][1] * m->v[3][2] -
                         m->v[2][0] * m->v[1][2] * m->v[3][1] -
                         m->v[3][0] * m->v[1][1] * m->v[2][2] +
                         m->v[3][0] * m->v[1][2] * m->v[2][1];

        inv_m->v[0][1] = -m->v[0][1] * m->v[2][2] * m->v[3][3] +
                         m->v[0][1] * m->v[2][3] * m->v[3][2] +
                         m->v[2][1] * m->v[0][2] * m->v[3][3] -
                         m->v[2][1] * m->v[0][3] * m->v[3][2] -
                         m->v[3][1] * m->v[0][2] * m->v[2][3] +
                         m->v[3][1] * m->v[0][3] * m->v[2][2];

        inv_m->v[1][1] = m->v[0][0] * m->v[2][2] * m->v[3][3] -
                         m->v[0][0] * m->v[2][3] * m->v[3][2] -
                         m->v[2][0] * m->v[0][2] * m->v[3][3] +
                         m->v[2][0] * m->v[0][3] * m->v[3][2] +
                         m->v[3][0] * m->v[0][2] * m->v[2][3] -
                         m->v[3][0] * m->v[0][3] * m->v[2][2];

        inv_m->v[2][1] = -m->v[0][0] * m->v[2][1] * m->v[3][3] +
                         m->v[0][0] * m->v[2][3] * m->v[3][1] +
                         m->v[2][0] * m->v[0][1] * m->v[3][3] -
                         m->v[2][0] * m->v[0][3] * m->v[3][1] -
                         m->v[3][0] * m->v[0][1] * m->v[2][3] +
                         m->v[3][0] * m->v[0][3] * m->v[2][1];

        inv_m->v[3][1] = m->v[0][0] * m->v[2][1] * m->v[3][2] -
                         m->v[0][0] * m->v[2][2] * m->v[3][1] -
                         m->v[2][0] * m->v[0][1] * m->v[3][2] +
                         m->v[2][0] * m->v[0][2] * m->v[3][1] +
                         m->v[3][0] * m->v[0][1] * m->v[2][2] -
                         m->v[3][0] * m->v[0][2] * m->v[2][1];

        inv_m->v[0][2] = m->v[0][1] * m->v[1][2] * m->v[3][3] -
                         m->v[0][1] * m->v[1][3] * m->v[3][2] -
                         m->v[1][1] * m->v[0][2] * m->v[3][3] +
                         m->v[1][1] * m->v[0][3] * m->v[3][2] +
                         m->v[3][1] * m->v[0][2] * m->v[1][3] -
                         m->v[3][1] * m->v[0][3] * m->v[1][2];

        inv_m->v[1][2] = -m->v[0][0] * m->v[1][2] * m->v[3][3] +
                         m->v[0][0] * m->v[1][3] * m->v[3][2] +
                         m->v[1][0] * m->v[0][2] * m->v[3][3] -
                         m->v[1][0] * m->v[0][3] * m->v[3][2] -
                         m->v[3][0] * m->v[0][2] * m->v[1][3] +
                         m->v[3][0] * m->v[0][3] * m->v[1][2];

        inv_m->v[2][2] = m->v[0][0] * m->v[1][1] * m->v[3][3] -
                         m->v[0][0] * m->v[1][3] * m->v[3][1] -
                         m->v[1][0] * m->v[0][1] * m->v[3][3] +
                         m->v[1][0] * m->v[0][3] * m->v[3][1] +
                         m->v[3][0] * m->v[0][1] * m->v[1][3] -
                         m->v[3][0] * m->v[0][3] * m->v[1][1];

        inv_m->v[3][2] = -m->v[0][0] * m->v[1][1] * m->v[3][2] +
                         m->v[0][0] * m->v[1][2] * m->v[3][1] +
                         m->v[1][0] * m->v[0][1] * m->v[3][2] -
                         m->v[1][0] * m->v[0][2] * m->v[3][1] -
                         m->v[3][0] * m->v[0][1] * m->v[1][2] +
                         m->v[3][0] * m->v[0][2] * m->v[1][1];

        inv_m->v[0][3] = -m->v[0][1] * m->v[1][2] * m->v[2][3] +
                         m->v[0][1] * m->v[1][3] * m->v[2][2] +
                         m->v[1][1] * m->v[0][2] * m->v[2][3] -
                         m->v[1][1] * m->v[0][3] * m->v[2][2] -
                         m->v[2][1] * m->v[0][2] * m->v[1][3] +
                         m->v[2][1] * m->v[0][3] * m->v[1][2];

        inv_m->v[1][3] = m->v[0][0] * m->v[1][2] * m->v[2][3] -
                         m->v[0][0] * m->v[1][3] * m->v[2][2] -
                         m->v[1][0] * m->v[0][2] * m->v[2][3] +
                         m->v[1][0] * m->v[0][3] * m->v[2][2] +
                         m->v[2][0] * m->v[0][2] * m->v[1][3] -
                         m->v[2][0] * m->v[0][3] * m->v[1][2];

        inv_m->v[2][3] = -m->v[0][0] * m->v[1][1] * m->v[2][3] +
                         m->v[0][0] * m->v[1][3] * m->v[2][1] +
                         m->v[1][0] * m->v[0][1] * m->v[2][3] -
                         m->v[1][0] * m->v[0][3] * m->v[2][1] -
                         m->v[2][0] * m->v[0][1] * m->v[1][3] +
                         m->v[2][0] * m->v[0][3] * m->v[1][1];

        inv_m->v[3][3] = m->v[0][0] * m->v[1][1] * m->v[2][2] -
                         m->v[0][0] * m->v[1][2] * m->v[2][1] -
                         m->v[1][0] * m->v[0][1] * m->v[2][2] +
                         m->v[1][0] * m->v[0][2] * m->v[2][1] +
                         m->v[2][0] * m->v[0][1] * m->v[1][2] -
                         m->v[2][0] * m->v[0][2] * m->v[1][1];

        double det = m->v[0][0] * inv_m->v[0][0] +
                     m->v[0][1] * inv_m->v[1][0] +
                     m->v[0][2] * inv_m->v[2][0] +
                     m->v[0][3] * inv_m->v[3][0];

        if (det == 0)
            return false;

        det = 1.0 / det;

        for (int i = 0; i < 4; i++)
        {
            for (int j = 0; j < 4; j++)
            {
                inv_m->v[i][j] = inv_m->v[i][j] * det;
            }
        }

        return true;
    }

Adding a 2d case which might be useful for someone else:

inline bool invert4by4matrix(mat m, mat inv_m)
 {

        inv_m->v[0][0] = m->v[1][1] * m->v[2][2] * m->v[3][3] -
                         m->v[1][1] * m->v[2][3] * m->v[3][2] -
                         m->v[2][1] * m->v[1][2] * m->v[3][3] +
                         m->v[2][1] * m->v[1][3] * m->v[3][2] +
                         m->v[3][1] * m->v[1][2] * m->v[2][3] -
                         m->v[3][1] * m->v[1][3] * m->v[2][2];

        inv_m->v[1][0] = -m->v[1][0] * m->v[2][2] * m->v[3][3] +
                         m->v[1][0] * m->v[2][3] * m->v[3][2] +
                         m->v[2][0] * m->v[1][2] * m->v[3][3] -
                         m->v[2][0] * m->v[1][3] * m->v[3][2] -
                         m->v[3][0] * m->v[1][2] * m->v[2][3] +
                         m->v[3][0] * m->v[1][3] * m->v[2][2];

        inv_m->v[2][0] = m->v[1][0] * m->v[2][1] * m->v[3][3] -
                         m->v[1][0] * m->v[2][3] * m->v[3][1] -
                         m->v[2][0] * m->v[1][1] * m->v[3][3] +
                         m->v[2][0] * m->v[1][3] * m->v[3][1] +
                         m->v[3][0] * m->v[1][1] * m->v[2][3] -
                         m->v[3][0] * m->v[1][3] * m->v[2][1];

        inv_m->v[3][0] = -m->v[1][0] * m->v[2][1] * m->v[3][2] +
                         m->v[1][0] * m->v[2][2] * m->v[3][1] +
                         m->v[2][0] * m->v[1][1] * m->v[3][2] -
                         m->v[2][0] * m->v[1][2] * m->v[3][1] -
                         m->v[3][0] * m->v[1][1] * m->v[2][2] +
                         m->v[3][0] * m->v[1][2] * m->v[2][1];

        inv_m->v[0][1] = -m->v[0][1] * m->v[2][2] * m->v[3][3] +
                         m->v[0][1] * m->v[2][3] * m->v[3][2] +
                         m->v[2][1] * m->v[0][2] * m->v[3][3] -
                         m->v[2][1] * m->v[0][3] * m->v[3][2] -
                         m->v[3][1] * m->v[0][2] * m->v[2][3] +
                         m->v[3][1] * m->v[0][3] * m->v[2][2];

        inv_m->v[1][1] = m->v[0][0] * m->v[2][2] * m->v[3][3] -
                         m->v[0][0] * m->v[2][3] * m->v[3][2] -
                         m->v[2][0] * m->v[0][2] * m->v[3][3] +
                         m->v[2][0] * m->v[0][3] * m->v[3][2] +
                         m->v[3][0] * m->v[0][2] * m->v[2][3] -
                         m->v[3][0] * m->v[0][3] * m->v[2][2];

        inv_m->v[2][1] = -m->v[0][0] * m->v[2][1] * m->v[3][3] +
                         m->v[0][0] * m->v[2][3] * m->v[3][1] +
                         m->v[2][0] * m->v[0][1] * m->v[3][3] -
                         m->v[2][0] * m->v[0][3] * m->v[3][1] -
                         m->v[3][0] * m->v[0][1] * m->v[2][3] +
                         m->v[3][0] * m->v[0][3] * m->v[2][1];

        inv_m->v[3][1] = m->v[0][0] * m->v[2][1] * m->v[3][2] -
                         m->v[0][0] * m->v[2][2] * m->v[3][1] -
                         m->v[2][0] * m->v[0][1] * m->v[3][2] +
                         m->v[2][0] * m->v[0][2] * m->v[3][1] +
                         m->v[3][0] * m->v[0][1] * m->v[2][2] -
                         m->v[3][0] * m->v[0][2] * m->v[2][1];

        inv_m->v[0][2] = m->v[0][1] * m->v[1][2] * m->v[3][3] -
                         m->v[0][1] * m->v[1][3] * m->v[3][2] -
                         m->v[1][1] * m->v[0][2] * m->v[3][3] +
                         m->v[1][1] * m->v[0][3] * m->v[3][2] +
                         m->v[3][1] * m->v[0][2] * m->v[1][3] -
                         m->v[3][1] * m->v[0][3] * m->v[1][2];

        inv_m->v[1][2] = -m->v[0][0] * m->v[1][2] * m->v[3][3] +
                         m->v[0][0] * m->v[1][3] * m->v[3][2] +
                         m->v[1][0] * m->v[0][2] * m->v[3][3] -
                         m->v[1][0] * m->v[0][3] * m->v[3][2] -
                         m->v[3][0] * m->v[0][2] * m->v[1][3] +
                         m->v[3][0] * m->v[0][3] * m->v[1][2];

        inv_m->v[2][2] = m->v[0][0] * m->v[1][1] * m->v[3][3] -
                         m->v[0][0] * m->v[1][3] * m->v[3][1] -
                         m->v[1][0] * m->v[0][1] * m->v[3][3] +
                         m->v[1][0] * m->v[0][3] * m->v[3][1] +
                         m->v[3][0] * m->v[0][1] * m->v[1][3] -
                         m->v[3][0] * m->v[0][3] * m->v[1][1];

        inv_m->v[3][2] = -m->v[0][0] * m->v[1][1] * m->v[3][2] +
                         m->v[0][0] * m->v[1][2] * m->v[3][1] +
                         m->v[1][0] * m->v[0][1] * m->v[3][2] -
                         m->v[1][0] * m->v[0][2] * m->v[3][1] -
                         m->v[3][0] * m->v[0][1] * m->v[1][2] +
                         m->v[3][0] * m->v[0][2] * m->v[1][1];

        inv_m->v[0][3] = -m->v[0][1] * m->v[1][2] * m->v[2][3] +
                         m->v[0][1] * m->v[1][3] * m->v[2][2] +
                         m->v[1][1] * m->v[0][2] * m->v[2][3] -
                         m->v[1][1] * m->v[0][3] * m->v[2][2] -
                         m->v[2][1] * m->v[0][2] * m->v[1][3] +
                         m->v[2][1] * m->v[0][3] * m->v[1][2];

        inv_m->v[1][3] = m->v[0][0] * m->v[1][2] * m->v[2][3] -
                         m->v[0][0] * m->v[1][3] * m->v[2][2] -
                         m->v[1][0] * m->v[0][2] * m->v[2][3] +
                         m->v[1][0] * m->v[0][3] * m->v[2][2] +
                         m->v[2][0] * m->v[0][2] * m->v[1][3] -
                         m->v[2][0] * m->v[0][3] * m->v[1][2];

        inv_m->v[2][3] = -m->v[0][0] * m->v[1][1] * m->v[2][3] +
                         m->v[0][0] * m->v[1][3] * m->v[2][1] +
                         m->v[1][0] * m->v[0][1] * m->v[2][3] -
                         m->v[1][0] * m->v[0][3] * m->v[2][1] -
                         m->v[2][0] * m->v[0][1] * m->v[1][3] +
                         m->v[2][0] * m->v[0][3] * m->v[1][1];

        inv_m->v[3][3] = m->v[0][0] * m->v[1][1] * m->v[2][2] -
                         m->v[0][0] * m->v[1][2] * m->v[2][1] -
                         m->v[1][0] * m->v[0][1] * m->v[2][2] +
                         m->v[1][0] * m->v[0][2] * m->v[2][1] +
                         m->v[2][0] * m->v[0][1] * m->v[1][2] -
                         m->v[2][0] * m->v[0][2] * m->v[1][1];

        double det = m->v[0][0] * inv_m->v[0][0] +
                     m->v[0][1] * inv_m->v[1][0] +
                     m->v[0][2] * inv_m->v[2][0] +
                     m->v[0][3] * inv_m->v[3][0];

        if (det == 0)
            return false;

        det = 1.0 / det;

        for (int i = 0; i < 4; i++)
        {
            for (int j = 0; j < 4; j++)
            {
                inv_m->v[i][j] = inv_m->v[i][j] * det;
            }
        }

        return true;
    }
回复 原文
孤君无依 2024-08-03 19:37:13

如果你想计算 4x4 矩阵的逆矩阵，那么我建议使用像 OpenGL数学（GLM）：

无论如何，你可以从头开始。 以下实现类似于 glm 的实现： :inverse，但它没有那么高度优化：

bool InverseMat44( const GLfloat m[16], GLfloat invOut[16] )
{
    float 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;
}

If you want to compute the inverse matrix of 4x4 matrix, then I recommend to use a library like OpenGL Mathematics (GLM) :

Anyway, you can do it from scratch. The following implementation is similar to the implementation of glm::inverse, but it is not as highly optimized:

bool InverseMat44( const GLfloat m[16], GLfloat invOut[16] )
{
    float 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;
}
回复 原文
~没有更多了~

关于作者

被你宠の有点坏

暂无简介

0 文章
0 评论
22 人气
发私信

相关话题

更多

热门标签

操作系统 程序设计 IT运维 Linux系统管理 JavaScript 服务器应用 solaris C/C++ PHP Shell BSD Vue.js aix Oracle Python HTML 系统管理 HTML5 CSS 前端
更多

推荐作者

内心激荡

文章 0 评论 0

JSmiles

文章 0 评论 0

赏烟花じ飞满天

文章 0 评论 0

左秋

文章 0 评论 0

迪街小绵羊

文章 0 评论 0

瞳孔里扚悲伤

文章 0 评论 0

更多

友情链接

文江博客
    我们使用 Cookies 和其他技术来定制您的体验包括您的登录状态等。通过阅读我们的 隐私政策 了解更多相关信息。 单击 接受 或继续使用网站，即表示您同意使用 Cookies 和您的相关数据。
    原文