đảo ngược ma trận 4x4

Question

tôi đang tìm cách triển khai mã mẫu về cách đảo ngược ma trận 4x4. tôi biết có gaussian eleminiation, LU phân hủy, vv nhưng thay vì nhìn vào chúng chi tiết tôi thực sự chỉ cần tìm kiếm mã để làm điều này.

ngôn ngữ lý tưởng C++, dữ liệu có sẵn trong mảng 16 phao theo thứ tự lớn.

cảm ơn bạn!

Answer 1

đây:

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

này đã được dỡ bỏ từ MESA triển khai thư viện GLU.

Answer 2

CHO 3x3 MATRIX

ĐỔI BỘ LUẬT THEO YÊU CẦU CỦA BẠN

http://www.dreamincode.net/code/snippet1156.htm

Cập nhật:

Có ... Đi từ 3x3 đến 4x4 có vẻ như một sự khác biệt lớn ... câu trả lời này là không chính xác cho việc này.

Answer 3

Bạn có thể sử dụng GNU Scientific Library hoặc tìm mã trong đó.

Chỉnh sửa: Bạn dường như muốn phần Linear Algebra.

Answer 4

Đây là thư viện nhỏ (chỉ một tiêu đề) C++ vector math (hướng đến lập trình 3D). Nếu bạn sử dụng nó, hãy nhớ rằng cách bố trí của các ma trận của nó trong bộ nhớ được đảo ngược so với những gì OpenGL hy vọng, tôi đã có thời gian vui vẻ tìm nó ra ...

Answer 5

Nếu bạn cần một Thư viện C++ ma trận với rất nhiều chức năng, có một cái nhìn tại Eigen thư viện - http://eigen.tuxfamily.org

Answer 6

tôi 'cuộn lại' việc thực hiện MESA (cũng viết một vài bài kiểm tra đơn vị để đảm bảo nó thực sự hoạt động).

đây:

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; 

} 

tôi đã viết một chút về vấn đề này và hiển thị các mô hình của/yếu tố tiêu cực tích cực on my blog.

Như đề xuất của @LiraNuna, trên nhiều nền tảng phần cứng tăng tốc các phiên bản của các thói quen như vậy có sẵn vì vậy tôi rất vui khi có một 'phiên bản dự phòng' có thể đọc được và súc tích.

Lưu ý: điều này có thể chạy chậm hơn hoặc chậm hơn 3,5 lần so với triển khai MESA. Bạn có thể thay đổi mô hình các yếu tố để loại bỏ một số bổ sung, v.v ... nhưng nó sẽ mất khả năng đọc và vẫn sẽ không nhanh.

Answer 7

Bạn có thể làm nhanh hơn theo blog này.

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

Đây không phải là một việc thực hiện hoàn chỉnh vì P có thể không khả nghịch, nhưng bạn có thể kết hợp với mã này thực hiện MESA để có được một hiệu suất tốt hơn.

Answer 8

Nếu bất cứ ai tìm mã costumized hơn và "dễ đọc", sau đó tôi nhận này

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 A= 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), 
}; 

Tôi không viết mã, nhưng chương trình của tôi đã làm. Tôi đã thực hiện một chương trình nhỏ để tạo một chương trình để tính toán yếu tố quyết định và nghịch đảo của bất kỳ ma trận N nào.

Tôi làm điều đó bởi vì một lần trong quá khứ tôi cần một mã mà đảo ngược ma trận 5x5, nhưng không ai trong trái đất đã làm điều này vì vậy tôi đã thực hiện một.

Hãy xem chương trình here.