高效的 4x4 矩阵逆（仿射变换）

Question

我希望有人能指出 4x4 仿射矩阵变换的有效公式。目前我的代码使用辅因子扩展，它为每个辅因子分配一个临时数组。它易于阅读，但速度比应有的要慢。

请注意，这不是家庭作业，我知道如何使用 4x4 协因子展开手动解决，这对我来说只是一个痛苦而不是真正有趣的问题。此外，我用谷歌搜索并想出了一些已经为您提供公式的网站（http://www.euclideanspace.com/maths/algebra/matrix/functions/inverse/fourD/index.htm）。然而，这可能可以通过预先计算一些产品来进一步优化。我敢肯定有人在某个时候为此想出了“最佳”公式？

Answer 1

我相信计算逆的唯一方法是解 n 次方程：A x = y，其中 y 跨越单位向量，即第一个是 (1,0,0,0)，第二个是（0,1,0,0）等

（使用辅因子（Cramer 规则）是个坏主意，除非你想要一个符号公式来表示逆。）

大多数线性代数库都可以让您求解这些线性系统，甚至可以计算逆。 python中的示例（使用numpy）：

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

Answer 2

IIRC 您可以通过预先计算一堆（12 个？）2x2 行列式来大大缩减代码和时间。将矩阵垂直分成两半，并在上半部分和下半部分每 2x2 计算一次。这些较小的决定因素之一用于您需要进行较大计算的每个术语，并且它们每个都可以重复使用。

另外，不要使用单独的行列式函数 - 重用您为伴随计算的子行列式来获得行列式。

哦，刚找到 this.

了解它的某种变换，您也可以做出一些改进。

Answer 3

您应该能够利用矩阵仿射这一事实来加速整个逆过程。也就是说，如果你的矩阵看起来像这样

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 向量（通常是 3D 点）。

最后，如果 M 表示旋转（即它的列是正交的），那么您可以使用 inv(M) = transpose(M) 的事实。那么计算 A 的逆只需减去平移分量，然后乘以 3x3 部分的转置即可。

请注意，矩阵是否正交是您应该从问题分析中知道的。在运行时检查它会相当昂贵；尽管您可能希望在调试版本中执行此操作以检查您的假设是否成立。

希望所有这些都清楚......

Answer 4

以防万一有人想节省一些打字，这是我根据上面 phkahler 发布的链接的第 9 页（拉普拉斯展开定理的更有效版本）编写的 AS3 版本：

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

当我将各种 3D 变换矩阵乘以从该方法返回的逆矩阵时，这成功生成了一个单位矩阵。我相信您可以搜索/替换以将其转换为您喜欢的任何语言。

Answer 5

为了跟进pkhaler 和Robin Hilliard 上面的出色回复，这里是Robin 的ActionScript 3 代码转换为C# 方法。希望这可以为其他 C# 开发人员以及需要 4x4 矩阵求逆函数的 C/C++ 和 Java 开发人员节省一些输入：

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