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

Question

Я надеялся, что кто-нибудь сможет указать эффективную формулу для аффинного матричного преобразования 4x4. В настоящее время мой код использует расширение кофактора, и он выделяет временный массив для каждого кофактора. Его легко читать, но он медленнее, чем должен быть.

Обратите внимание, это не домашнее задание, и я знаю, как это сделать вручную, используя расширение с коэффициентом 4x4, это просто боль и не очень интересная проблема для меня. Кроме того, я гуглил и придумал несколько сайтов, которые уже дают формулу (http://www.euclideanspace.com/maths/algebra/matrix/functions/inverse/fourD/index.htm). Однако, возможно, этот вариант можно еще больше оптимизировать, предварительно рассчитав некоторые продукты. «лучшая» формула для этого в тот или иной момент?

Answer 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). Тогда вычисление обратной величины A - это всего лишь вопрос вычитания компонента перевода и умножения на транспонирование части 3x3.

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

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

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

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

Answer 3

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

Answer 4

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

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

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

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

Answer 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