Efficace 4x4 matrice inverse (affine transform)

j'espérais que quelqu'un puisse faire remarquer une formule efficace pour la transformation de matrice 4x4 affine. Actuellement, my code utilise l'extension du cofacteur et attribue un tableau temporaire pour chaque cofacteur. C'est facile à lire, mais c'est plus lent que prévu.

Note, Ce n'est pas un devoir et je sais comment le faire manuellement en utilisant l'expansion 4x4 co-facteur, c'est juste une douleur et pas vraiment un problème intéressant pour moi. J'ai aussi cherché sur Google et j'ai trouvé quelques sites qui vous donnent le formule déjà (http://www.euclideanspace.com/maths/algebra/matrix/functions/inverse/fourD/index.htm). Toutefois, celui-ci pourrait probablement être optimisé davantage en pré-calculant certains des produits. Je suis sûr que quelqu'un a trouvé la "meilleure" formule pour ça à un moment ou à un autre?

26
demandé sur vaxquis 2010-04-12 22:33:50

5 réponses

Vous devriez être en mesure d'exploiter le fait que la matrice est affine pour accélérer les choses pour l'ensemble de l'inverse. À savoir, si votre matrice ressemble à ceci

A = [ M   b  ]
    [ 0   1  ]

où Un 4x4, M 3x3, b est 3x1, et la rangée du bas (0,0,0,1), puis

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

selon votre situation, il peut être plus rapide de calculer le résultat de inv(A) * x au lieu de former réellement inv(A). Dans ce cas, les choses se simplifient

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

où x est un vecteur 3x1 (généralement un vecteur 3D point.)

enfin, si M représente une rotation(ses colonnes sont orthonormales), alors vous pouvez utiliser le fait que inv(M) = transpose (M). Ensuite, calculer L'inverse de A n'est qu'une question de soustraire la composante de traduction, et de multiplier par la transposition de la partie 3x3.

notez que si oui ou non la matrice est orthonormale est quelque chose que vous devriez savoir de l'analyse du problème. Le vérifier pendant l'exécution serait assez coûteux; bien que vous pourriez voulez le faire en debug pour vérifier vos hypothèses.

j'Espère tout cela est clair...

41
répondu celion 2010-04-13 18:19:51

juste au cas où quelqu'un voudrait sauver un peu de Dactylographie, voici une version AS3 que j'ai écrit basé sur la page 9 (version plus efficace du théorème D'Expansion de Laplace) du lien posté ci-dessus par 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;
}

cela a produit avec succès une matrice d'identité lorsque j'ai multiplié diverses matrices de transformation 3D par l'inverse retourné de cette méthode. Je suis sûr que vous pouvez rechercher/remplacer pour obtenir ce dans quelle langue vous le souhaitez.

19
répondu Robin Hilliard 2013-06-14 22:27:04

suivi pkhaleret Robin Hilliard's excellentes réponses ci-dessus, voici le code ActionScript 3 de Robin converti en une méthode C#. Avec un peu de chance, cela permettra de sauvegarder des caractères pour d'autres développeurs C#, ainsi que pour les développeurs C/C++ et Java qui ont besoin d'une fonction d'inversion de 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;
}
16
répondu Anders Gustafsson 2017-05-23 10:30:43

IIRC vous pouvez réduire considérablement le code et le temps en précomptant un tas (12?) 2x2 déterminants. Diviser la matrice en deux verticalement et calculer tous les 2x2 dans la moitié supérieure et inférieure. Un de ces déterminants plus petits est utilisé dans chaque terme que vous aurez besoin pour le calcul plus grand et ils sont chacun réutilisés.

aussi, n'utilisez pas une fonction de déterminant séparée - réutilisez les sous-déterminants que vous avez calculés pour l'adjoint pour obtenir le déterminant.

Oh, juste trouvé ceci.

il y a quelques améliorations que vous pouvez faire savoir que c'est aussi un certain type de transformation.

4
répondu phkahler 2010-04-12 18:56:27

je crois que la seule façon de calculer un inverse est de résoudre n fois l'équation: A x = y, où y couvre les vecteurs unitaires, i.e., le premier est (1.0.0.0), le second est (0.1.0.0), etc.

(utiliser les cofacteurs (règle de Cramer) est une mauvaise idée, sauf si vous voulez une formule symbolique pour l'inverse.)

la plupart des bibliothèques d'algèbre linéaire vous permettra de résoudre ces systèmes linéaires, et même de calculer un inverse. Exemple en python (utilisant numpy):

from numpy.linalg import inv
inv(A) # here you go
1
répondu Olivier Verdier 2010-04-12 18:55:41