Эффективная матрица 4x4 обратная (аффинное преобразование)

Ответ 1

Вы должны иметь возможность использовать тот факт, что матрица является аффинной, чтобы ускорить работу над полным обратным. А именно, если ваша матрица выглядит так:

A = [ M   b  ]
    [ 0   1  ]

где A - 4x4, M - 3x3, b - 3x1, а нижняя строка - (0,0,0,1), затем

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

В зависимости от вашей ситуации, скорее всего, быстрее вычислять результат inv (A) * x вместо фактического формирования inv (A). В этом случае вещи упрощаются до

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

где x - вектор 3x1 (обычно трехмерная точка).

Наконец, если M представляет собой поворот (т.е. его столбцы являются ортонормированными), то вы можете использовать тот факт, что inv (M) = transpose (M). Тогда вычисление обратного к А является просто вычитанием компонента трансляции и умножением на транспонирование части 3х3.

Обратите внимание, что независимо от того, является ли эта матрица ортонормированной, это то, что вы должны знать из анализа проблемы. Проверка его во время выполнения была бы довольно дорогой; хотя вы можете захотеть сделать это в отладочных сборках, чтобы проверить, что ваши предположения сохраняются.

Надеюсь, все это ясно...

Ответ 2

На всякий случай кто-то хотел бы сохранить некоторую типизацию, вот версия AS3, которую я написал на основе страницы 9 (более эффективная версия теоремы расширения Лапласа) ссылки, размещенной выше 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;
}

Это успешно создало единичную матрицу, когда я умножил различные матрицы трехмерного преобразования на обратный, возвращенный из этого метода. Я уверен, что вы можете искать/заменять, чтобы получить это на любом языке, который вам нужен.

Ответ 3

Чтобы следить за pkhaler и Робин Хиллиард, отличные ответы выше, вот код Robin ActionScript 3 преобразованный в метод С#. Надеюсь, это может сохранить некоторые типизации для других разработчиков С#, а также разработчиков C/С++ и Java, нуждающихся в функции инверсии матрицы 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;
}

Ответ 4

IIRC вы можете значительно сократить код и время, предварительно вычислив детерминанты 2x2 (12?) 2x2. Разделите матрицу пополам по вертикали и вычислите каждые 2x2 как в верхней, так и в нижней половине. Один из этих меньших детерминант используется во всех терминах, которые вам понадобятся для больших вычислений, и каждый из них будет использоваться повторно.

Кроме того, не используйте отдельную функцию детерминанта - повторно используйте подспределители, которые вы вычислили для сопряженного, чтобы получить определитель.

О, только что нашел this.

Есть некоторые улучшения, которые вы можете сделать, зная его определенный вид преобразования.

Ответ 5

Я считаю, что единственный способ вычислить обратный - это решить n раз уравнение: A x = y, где y охватывает единичные векторы, т.е. первый из них (1,0,0,0), второй (0,1,0,0) и т.д.

(Использование кофакторов (правило Крамера) - плохая идея, если вы не хотите использовать символическую формулу для обратного.)

Большинство библиотек линейной алгебры позволят вам решать эти линейные системы и даже вычислять обратные. Пример в python (с использованием numpy):

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