Efektywna odwrotność macierzy 4x4 (transformacja afiniczna)

Question

Miałem nadzieję, że ktoś może wskazać skuteczną formułę transformacji macierzy afinicznej 4x4. Obecnie mój kod wykorzystuje rozszerzenie kofaktora i przydziela tymczasową tablicę dla każdego kofaktora. Jest łatwy do odczytania, ale wolniejszy niż powinien.

Uwaga, to nie jest praca domowa i wiem, jak to zrobić ręcznie za pomocą rozszerzenia kolidującego 4x4, to tylko ból i nie jest to naprawdę interesujący problem dla mnie. Również googlowałem i wymyśliłem kilka stron, które podają już wzór (http://www.euclideanspace.com/maths/algebra/matrix/functions/inverse/fourD/index.htm). Jednak ta prawdopodobnie mogłaby zostać zoptymalizowana poprzez wstępne obliczenia niektórych produktów. Jestem pewien, że ktoś wymyślił "najlepszą" formułę do tego w takim czy innym punkcie?

Answer 1

Powinieneś być w stanie wykorzystać fakt, że matryca jest podatna na przyspieszenie rzeczy na pełną odwrotność. Mianowicie, jeśli matryca wygląda tak

A = [ M b ] 
    [ 0 1 ] 

gdzie A jest 4x4, 3x3 M, B jest 3x1, a dolny rząd jest (0,0,0,1), a następnie

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

W zależności od sytuacji może być szybsze wyliczenie wyniku inv (A) * x zamiast faktycznego utworzenia inv (A). W takim przypadku rzeczy upraszczają do

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

gdzie x to wektor 3x1 (zwykle punkt 3D).

Wreszcie, jeśli M reprezentuje obrót (to znaczy, że jego kolumny są ortonormalne), wówczas można użyć faktu, że inv (M) = transponuj (M). Następnie obliczenie odwrotności A jest tylko kwestią odjęcia komponentu translacji i pomnożenie przez transpozycję części 3x3.

Należy zauważyć, że to, czy macierz jest ortonormalna, jest czymś, o czym powinieneś wiedzieć z analizy problemu. Sprawdzenie go w czasie działania byłoby dość kosztowne; chociaż możesz chcieć to zrobić w kompilacjach debugujących, aby sprawdzić, czy twoje założenia są spełnione.

Nadzieja to wszystko jest jasne ...

Answer 2

Wierzę, że jedynym sposobem obliczenia odwrotności jest rozwiązanie n razy równania: A x = y, gdzie y obejmuje wektory jednostek, tj. Pierwszy to (1,0,0,0), drugi (0,1,0,0), itp

(za pomocą kofaktorów (Wzory Cramera) jest zły pomysł, chyba że chcesz symboliczny wzór na odwrotność.)

Większość bibliotek algebry liniowej pozwoli ci rozwiązać te systemy liniowe, a nawet obliczyć odwrotność. Przykład w Pythonie (przy użyciu numpy):

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

Answer 3

IIRC można znacznie zmniejszyć kod i czas precomputing pęczek (12?) Determinanty 2x2. Podziel macierz na pół w pionie i obliczaj co 2x2 w górnej i dolnej połowie. Jedna z tych mniejszych determinant jest używana w każdym terminie potrzebnym do większych obliczeń i każdy z nich jest ponownie wykorzystywany.

Ponadto nie należy używać oddzielnej funkcji wyznaczania - należy ponownie wykorzystać wyznaczone pod kątem determinanty obliczenia dla łącznika, aby uzyskać wyznacznik.

Och, po prostu znaleźć this.

Istnieją pewne ulepszenia można zrobić znając jej pewnego rodzaju przekształcać też.

Answer 4

Tylko w przypadku, gdy ktoś chciałby zaoszczędzić trochę pisać, oto wersja AS3 pisałem na podstawie strony 9 (bardziej wydajna wersja twierdzenia Expansion Laplace'a) z linkiem wysłane powyżej 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; 
} 

ten pomyślnie wytwarza macierz tożsamości, gdy pomnożony różne 3D macierzy transformacji przez odwrotną zwrócony z tego sposobu. Jestem pewien, że możesz wyszukiwać/zastępować, aby uzyskać to w dowolnym języku.

Answer 5

W celu sprawdzenia doskonałych odpowiedzi powyżej pkhaler i Robin Hilliard, tutaj kod ActionScript 3 Robin'a jest konwertowany na metodę C#. Miejmy nadzieję, że pozwoli to zaoszczędzić trochę pisania dla innych programistów C#, a także dla programistów C/C++ i Java potrzebujących funkcji inwersji macierzy 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; 
}