Войти
ПрограммированиеПодсказкиОбщее

Инверсия матрицы 4x4

Автор:

Ниже приведен вариант инверсии для матриц 4x4, где порядок элементов строки горизонтальный (RowMajor). Не самый быстрый, но вполне рабочий (проверено под OpenGL) код для Pascal/Delphi.

type
  PMatrix3 = ^TMatrix3;
  PMatrix4 = ^TMatrix4;
  TVector3 = array[0..2] of Single;
  TVector4 = array[0..3] of Single;
  TMatrix3 = array[0..2] of TVector3;
  TMatrix4 = array[0..3] of TVector4;

function GetMinor4(var M: TMatrix4; ACol, ARow: Integer): Single;
var
  R, C: Integer;
  N: TMatrix3;
  NC, NR: Integer;

begin
  NR := 0; // Строки
  NC := 0; // Столбцы

  // Строки
  for R := 0 to 3 do
  begin
    // Исключаем строку матрицы
    if (R <> ARow) then
    begin
      // Столбцы
      for C := 0 to 3 do
      begin
        // Исключаем столбец матрицы
        if (C <> ACol) then
        begin
          N[NC, NR] := M[C, R];
          NR := (NR + 1);
        end;
      end;
      NC := (NC + 1);
      NR := 0;
    end;//if
  end;//for C

  // Result = a0*b1*c2 + a2*b0*c1 + a1*b2*c0 - a2*b1*c0 - a0*b2*c1 - a1*b0*c2
  Result := N[0, 0] * N[1, 1] * N[2, 2] +
            N[2, 0] * N[0, 1] * N[1, 2] +
            N[1, 0] * N[2, 1] * N[0, 2] -
            N[2, 0] * N[1, 1] * N[0, 2] -
            N[0, 0] * N[2, 1] * N[1, 2] -
            N[1, 0] * N[0, 1] * N[2, 2];
end;

function InverseMatrix4(M: TMatrix4): TMatrix4;
var
  detA, invA: Single;
  RM: TMatrix4 absolute Result;
  c, r: Integer;

begin
  // 1. Находим определитель матрицы 4x4
  detA := M[0, 0] * GetMinor4(M, 0, 0) -
          M[0, 1] * GetMinor4(M, 0, 1) +
          M[0, 2] * GetMinor4(M, 0, 2) -
          M[0, 3] * GetMinor4(M, 0, 3);

  // Обратная матрица существует при (detA ≠ 0)
  if (detA <> 0) then
  begin
    // Общий множитель для матрицы
    invA := (1.0 / detA);

    // 2. Находим матрицу миноров
    // 3. и учитываем матрицу алгебраических дополнений
    // [+, -, +, -]
    // [-, +, -, +]
    // [+, -, +, -]
    // [-, +, -, +]
    RM[0, 0] := GetMinor4(M, 0, 0) * invA;
    RM[1, 0] := -GetMinor4(M, 0, 1) * invA;
    RM[2, 0] := GetMinor4(M, 0, 2) * invA;
    RM[3, 0] := -GetMinor4(M, 0, 3) * invA;
    RM[0, 1] := -GetMinor4(M, 1, 0) * invA;
    RM[1, 1] := GetMinor4(M, 1, 1) * invA;
    RM[2, 1] := -GetMinor4(M, 1, 2) * invA;
    RM[3, 1] := GetMinor4(M, 1, 3) * invA;
    RM[0, 2] := GetMinor4(M, 2, 0) * invA;
    RM[1, 2] := -GetMinor4(M, 2, 1) * invA;
    RM[2, 2] := GetMinor4(M, 2, 2) * invA;
    RM[3, 2] := -GetMinor4(M, 2, 3) * invA;
    RM[0, 3] := -GetMinor4(M, 3, 0) * invA;
    RM[1, 3] := GetMinor4(M, 3, 1) * invA;
    RM[2, 3] := -GetMinor4(M, 3, 2) * invA;
    RM[3, 3] := GetMinor4(M, 3, 3) * invA;
  end
  else RM := M;
end;

Скорость для Intel Core i7-6950X 3 ГГц:
Pascal-версия — 1 167 024 вызовов в секунду (~2570 тактов).
SIMD-версия — 3 716 643 вызовов в секунду (~800 тактов).

Возможно существуют более быстрые способы инверсии,
но на данный момент мне они неизвестны.

Всем удачи!

25 января 2020 (Обновление: 28 янв. 2020)

Комментарии [86]