Efficiente matrice 4x4 inversa (trasformazione affine)

Question

Speravo che qualcuno potesse indicare una formula efficiente per la trasformazione della matrice affine 4x4. Attualmente il mio codice utilizza l'espansione cofattore e assegna un array temporaneo per ogni cofattore. È facile da leggere, ma è più lento di quanto dovrebbe essere.

Nota, questo non è compito a casa e so come lavorarlo manualmente utilizzando l'espansione co-factor 4x4, è solo un problema per me non particolarmente interessante. Anche io ho cercato su Google e ho trovato alcuni siti che ti danno già la formula (http://www.euclideanspace.com/maths/algebra/matrix/functions/inverse/fourD/index.htm). Tuttavia, questo potrebbe probabilmente essere ulteriormente ottimizzato pre-calcolando alcuni dei prodotti. Sono sicuro che qualcuno ha inventato la formula "migliore" per questo in un punto o nell'altro?

Answer 1

Si dovrebbe essere in grado di sfruttare il fatto che la matrice è affine per accelerare le cose su un inverso completo. Vale a dire, se la matrice si presenta così

A = [ M b ] 
    [ 0 1 ] 

dove A è 4x4, 3x3 M è, b è 3x1, e la fila inferiore è (0,0,0,1), quindi

inv(A) = [ inv(M) -inv(M) * b ] 
     [ 0   1  ] 

A seconda della situazione, potrebbe essere più veloce calcolare il risultato di inv (A) * x invece di formare effettivamente inv (A). In tal caso, le cose si semplificano in

inv(A) * [x] = [ inv(M) * (x - b) ] 
     [1] = [  1   ] 

dove x è un vettore 3x1 (di solito un punto 3D).

Infine, se M rappresenta una rotazione (cioè le sue colonne sono ortonormali), è possibile utilizzare il fatto che inv (M) = traspone (M). Quindi calcolare l'inverso di A è solo questione di sottrarre il componente di traduzione e moltiplicare per la trasposizione della parte 3x3.

Si noti che la matrice è ortonormale o qualcosa che dovresti sapere dall'analisi del problema. Controllarlo durante il runtime sarebbe abbastanza costoso; anche se potresti volerlo fare in build di debug per verificare che le tue ipotesi siano valide.

Speranza tutto questo è chiaro ...

Answer 2

Credo che l'unico modo per calcolare un inverso è quello di risolvere n volte l'equazione: A x = y, dove si estende i vettori unitari, cioè, il primo è (1,0,0,0), il in secondo luogo è (0,1,0,0), ecc

(utilizzando i cofattori (regola di Cramer) è una cattiva idea, a meno che non si desidera una formula simbolica per l'inverso.)

maggior parte delle librerie di algebra lineare ti permetterà di risolvere quei sistemi lineari e persino di calcolare un inverso. Esempio in Python (usando NumPy):

from numpy.linalg import inv 
inv(A) # here you go 

Answer 3

IIRC si può notevolmente ridurre il codice e tempo precomputing un gruppo (12?) Determinanti 2x2. Dividi la matrice a metà verticalmente e calcola ogni 2x2 sia nella metà superiore che in quella inferiore. Uno di questi determinanti minori viene utilizzato in tutti i termini necessari per il calcolo più grande e ognuno di essi viene riutilizzato.

Inoltre, non utilizzare una funzione determinante separata: riutilizzare i determinanti secondari calcolati per l'aggiornamento per ottenere il determinante.

Oh, appena trovato this.

ci sono alcuni miglioramenti si possono fare conoscere la sua un certo tipo di trasformare anche.

Answer 4

Solo nel caso qualcuno vorrebbe risparmiare un po 'di battitura, ecco una versione AS3 ho scritto in base a pagina 9 (versione più efficiente di Laplace espansione Teorema) del collegamento postato sopra per phkahler:

public function invert() : Matrix4 { 
    var m : Matrix4 = new Matrix4(); 

    var s0 : Number = i00 * i11 - i10 * i01; 
    var s1 : Number = i00 * i12 - i10 * i02; 
    var s2 : Number = i00 * i13 - i10 * i03; 
    var s3 : Number = i01 * i12 - i11 * i02; 
    var s4 : Number = i01 * i13 - i11 * i03; 
    var s5 : Number = i02 * i13 - i12 * i03; 

    var c5 : Number = i22 * i33 - i32 * i23; 
    var c4 : Number = i21 * i33 - i31 * i23; 
    var c3 : Number = i21 * i32 - i31 * i22; 
    var c2 : Number = i20 * i33 - i30 * i23; 
    var c1 : Number = i20 * i32 - i30 * i22; 
    var c0 : Number = i20 * i31 - i30 * i21; 

    // Should check for 0 determinant 

    var invdet : Number = 1/(s0 * c5 - s1 * c4 + s2 * c3 + s3 * c2 - s4 * c1 + s5 * c0); 

    m.i00 = (i11 * c5 - i12 * c4 + i13 * c3) * invdet; 
    m.i01 = (-i01 * c5 + i02 * c4 - i03 * c3) * invdet; 
    m.i02 = (i31 * s5 - i32 * s4 + i33 * s3) * invdet; 
    m.i03 = (-i21 * s5 + i22 * s4 - i23 * s3) * invdet; 

    m.i10 = (-i10 * c5 + i12 * c2 - i13 * c1) * invdet; 
    m.i11 = (i00 * c5 - i02 * c2 + i03 * c1) * invdet; 
    m.i12 = (-i30 * s5 + i32 * s2 - i33 * s1) * invdet; 
    m.i13 = (i20 * s5 - i22 * s2 + i23 * s1) * invdet; 

    m.i20 = (i10 * c4 - i11 * c2 + i13 * c0) * invdet; 
    m.i21 = (-i00 * c4 + i01 * c2 - i03 * c0) * invdet; 
    m.i22 = (i30 * s4 - i31 * s2 + i33 * s0) * invdet; 
    m.i23 = (-i20 * s4 + i21 * s2 - i23 * s0) * invdet; 

    m.i30 = (-i10 * c3 + i11 * c1 - i12 * c0) * invdet; 
    m.i31 = (i00 * c3 - i01 * c1 + i02 * c0) * invdet; 
    m.i32 = (-i30 * s3 + i31 * s1 - i32 * s0) * invdet; 
    m.i33 = (i20 * s3 - i21 * s1 + i22 * s0) * invdet; 

    return m; 
} 

questo prodotto con successo una matrice di identità quando ho moltiplicato varie matrici di trasformazione 3D per l'inverso restituito da questo metodo. Sono sicuro che puoi cercare/sostituire per ottenere ciò in qualsiasi lingua desideri.

Answer 5

Per dare seguito alle eccellenti risposte di pkhaler e qui sopra, ecco il codice ActionScript 3 di Robin convertito in un metodo C#. Speriamo che questo può risparmiare un po 'di battitura per gli altri sviluppatori C#, così come C/C++ e Java gli sviluppatori hanno bisogno di una funzione di inversione di matrice 4x4:

public static double[,] GetInverse(double[,] a) 
{ 
    var s0 = a[0, 0] * a[1, 1] - a[1, 0] * a[0, 1]; 
    var s1 = a[0, 0] * a[1, 2] - a[1, 0] * a[0, 2]; 
    var s2 = a[0, 0] * a[1, 3] - a[1, 0] * a[0, 3]; 
    var s3 = a[0, 1] * a[1, 2] - a[1, 1] * a[0, 2]; 
    var s4 = a[0, 1] * a[1, 3] - a[1, 1] * a[0, 3]; 
    var s5 = a[0, 2] * a[1, 3] - a[1, 2] * a[0, 3]; 

    var c5 = a[2, 2] * a[3, 3] - a[3, 2] * a[2, 3]; 
    var c4 = a[2, 1] * a[3, 3] - a[3, 1] * a[2, 3]; 
    var c3 = a[2, 1] * a[3, 2] - a[3, 1] * a[2, 2]; 
    var c2 = a[2, 0] * a[3, 3] - a[3, 0] * a[2, 3]; 
    var c1 = a[2, 0] * a[3, 2] - a[3, 0] * a[2, 2]; 
    var c0 = a[2, 0] * a[3, 1] - a[3, 0] * a[2, 1]; 

    // Should check for 0 determinant 
    var invdet = 1.0/(s0 * c5 - s1 * c4 + s2 * c3 + s3 * c2 - s4 * c1 + s5 * c0); 

    var b = new double[4, 4]; 

    b[0, 0] = (a[1, 1] * c5 - a[1, 2] * c4 + a[1, 3] * c3) * invdet; 
    b[0, 1] = (-a[0, 1] * c5 + a[0, 2] * c4 - a[0, 3] * c3) * invdet; 
    b[0, 2] = (a[3, 1] * s5 - a[3, 2] * s4 + a[3, 3] * s3) * invdet; 
    b[0, 3] = (-a[2, 1] * s5 + a[2, 2] * s4 - a[2, 3] * s3) * invdet; 

    b[1, 0] = (-a[1, 0] * c5 + a[1, 2] * c2 - a[1, 3] * c1) * invdet; 
    b[1, 1] = (a[0, 0] * c5 - a[0, 2] * c2 + a[0, 3] * c1) * invdet; 
    b[1, 2] = (-a[3, 0] * s5 + a[3, 2] * s2 - a[3, 3] * s1) * invdet; 
    b[1, 3] = (a[2, 0] * s5 - a[2, 2] * s2 + a[2, 3] * s1) * invdet; 

    b[2, 0] = (a[1, 0] * c4 - a[1, 1] * c2 + a[1, 3] * c0) * invdet; 
    b[2, 1] = (-a[0, 0] * c4 + a[0, 1] * c2 - a[0, 3] * c0) * invdet; 
    b[2, 2] = (a[3, 0] * s4 - a[3, 1] * s2 + a[3, 3] * s0) * invdet; 
    b[2, 3] = (-a[2, 0] * s4 + a[2, 1] * s2 - a[2, 3] * s0) * invdet; 

    b[3, 0] = (-a[1, 0] * c3 + a[1, 1] * c1 - a[1, 2] * c0) * invdet; 
    b[3, 1] = (a[0, 0] * c3 - a[0, 1] * c1 + a[0, 2] * c0) * invdet; 
    b[3, 2] = (-a[3, 0] * s3 + a[3, 1] * s1 - a[3, 2] * s0) * invdet; 
    b[3, 3] = (a[2, 0] * s3 - a[2, 1] * s1 + a[2, 2] * s0) * invdet; 

    return b; 
}