反转一个4x4matrix
我正在寻找如何反转4x4matrix的示例代码实现。 我知道有高斯eleminiation，LU分解等，而不是看他们的细节，我真的只是寻找代码来做到这一点。
语言是理想的C ++，数据以16个浮点数的forms以主要的顺序可用。
谢谢！
这里：
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; } 这是从MESA实施的GLU图书馆中解除的。
如果您正在寻找也是非常优化的“正常工作”实现，而不需要理解代码，我强烈推荐使用这里描述的英特尔优化的SSEmatrix逆向程序 。 在C中还有一个参考强制高斯消除和克莱默的规则。
我警告，如果你不了解MMX / SSE内在函数， 那么SSE代码就不好看了。
如果你需要一个有很多函数的C ++matrix库，看看Eigen库 – http://eigen.tuxfamily.org
我“汇总”了MESA实现（也编写了几个unit testing来确保实际工作）。
这里：
float invf(int i,int j,const float* m){ int o = 2+(ji); 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所build议的那样，在许多平台上，这些例程的硬件加速版本是可用的，所以我很高兴有一个可读和简洁的“备份版本”。
注意 ：这可能比MESA执行慢3.5倍或更差。 你可以改变因素的模式，以删除一些增加等等，但它会失去可读性，但仍然不会很快。
这是一个小的（只有一个头）C ++ vectormath库（面向3D编程）。 如果你使用它，记住它的matrix在内存中的布局与OpenGL期望的相反，我玩得很开心。
您可以使用GNU科学图书馆或查看代码。
编辑：你似乎想要线性代数部分。
根据这个博客，你可以加快速度。
#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实现结合起来以获得更好的性能。
如果有人正在寻找更多的代码和“更容易阅读”，那么我得到了这个
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 ), };
我不写代码，但我的程序。 我做了一个小程序来制作一个计算任何Nmatrix的行列式和逆matrix的程序。
我这样做是因为一旦过去我需要一个反转5x5matrix的代码，但地球上没有人做到这一点，所以我做了一个。
看看这里的程序。
对于3x3matrix
根据您的要求更改代码
http://www.dreamincode.net/code/snippet1156.htm
更新：
是的…从3×3到4×4似乎是一个很大的区别…这个答案是不正确的。
