-
Все привет! Кто-нибудь знает как инвертировать матрицу 4x3?
3x3 сделал, а 4x3 никак. У меня по методу Крамера. Формул для 4x4 ни у кого не завалялось?
-
Метод Крамера работает только с однозначно определёнными матрицами.
А матрица 4х3 гипотетически переопределённая. Для неё в общем случае считается псевдо-обртаная матрица. Приближёнными методами к примеру через SVD. А вообще возьмите нормальный движок и посмотрите как там сделано. Нормальный это OGRE и id Tech 3.
4х3 не обращают. Берут афинну часть 3х3 её обращают и знак у смещения меняют. Правда это требует поддержания афинности матрицы, но для игр и 3 графики это несложно сделать.
Вроде у меня обращение через миноры для 4х4
function Invert(A:TMatrix44):TMatrix44;Overload;
var det:real;
begin
det:=Determinant(A);
if det=0 then
begin
Result:=ZeroMatrix44;
end else
begin
det:=1/det;
Result[0,0]:=+det*(A[1,1]*(A[2,2]*A[3,3]-A[3,2]*A[2,3])-
A[1,2]*(A[2,1]*A[3,3]-A[3,1]*A[2,3])+
A[1,3]*(A[2,1]*A[3,2]-A[3,1]*A[2,2]));
Result[1,0]:=-det*(A[1,0]*(A[2,2]*A[3,3]-A[3,2]*A[2,3])-
A[1,2]*(A[2,0]*A[3,3]-A[3,0]*A[2,3])+
A[1,3]*(A[2,0]*A[3,2]-A[3,0]*A[2,2]));
Result[2,0]:=+det*(A[1,0]*(A[2,1]*A[3,3]-A[3,1]*A[2,3])-
A[1,1]*(A[2,0]*A[3,3]-A[3,0]*A[2,3])+
A[1,3]*(A[2,0]*A[3,1]-A[3,0]*A[2,1]));
Result[3,0]:=-det*(A[1,0]*(A[2,1]*A[3,2]-A[3,1]*A[2,2])-
A[1,1]*(A[2,0]*A[3,2]-A[3,0]*A[2,2])+
A[1,2]*(A[2,0]*A[3,1]-A[3,0]*A[2,1]));
Result[0,1]:=-det*(A[0,1]*(A[2,2]*A[3,3]-A[3,1]*A[2,3])-
A[0,2]*(A[2,1]*A[3,3]-A[3,1]*A[2,3])+
A[0,3]*(A[2,1]*A[3,2]-A[3,1]*A[2,2]));
Result[0,2]:=+det*(A[0,1]*(A[1,2]*A[3,3]-A[3,2]*A[1,3])-
A[0,2]*(A[1,1]*A[3,3]-A[3,1]*A[1,3])+
A[0,3]*(A[1,1]*A[3,2]-A[3,1]*A[1,2]));
Result[0,3]:=-det*(A[0,1]*(A[1,2]*A[2,3]-A[2,2]*A[1,3])-
A[0,2]*(A[1,1]*A[2,3]-A[2,1]*A[1,3])+
A[0,3]*(A[1,1]*A[2,2]-A[2,1]*A[1,2]));
Result[1,1]:=+det*(A[0,0]*(A[2,2]*A[3,3]-A[3,2]*A[2,3])-
A[0,2]*(A[2,0]*A[3,3]-A[3,0]*A[2,3])+
A[0,3]*(A[2,0]*A[3,2]-A[3,0]*A[2,2]));
Result[1,2]:=-det*(A[0,0]*(A[1,2]*A[3,3]-A[3,2]*A[1,3])-
A[0,2]*(A[1,0]*A[3,3]-A[3,0]*A[1,3])+
A[0,3]*(A[1,0]*A[3,2]-A[3,0]*A[1,2]));
Result[1,3]:=+det*(A[0,0]*(A[1,2]*A[2,3]-A[2,2]*A[1,3])-
A[0,2]*(A[1,0]*A[2,3]-A[2,0]*A[1,3])+
A[0,3]*(A[1,0]*A[2,2]-A[2,0]*A[1,2]));
Result[2,1]:=-det*(A[0,0]*(A[2,1]*A[3,3]-A[3,1]*A[2,3])-
A[0,1]*(A[2,0]*A[3,3]-A[3,0]*A[2,3])+
A[0,3]*(A[2,0]*A[3,1]-A[3,0]*A[2,1]));
Result[2,2]:=+det*(A[0,0]*(A[1,1]*A[3,3]-A[3,1]*A[1,3])-
A[0,1]*(A[1,0]*A[3,3]-A[3,0]*A[1,3])+
A[0,3]*(A[1,0]*A[3,1]-A[3,0]*A[1,1]));
Result[2,3]:=-det*(A[0,0]*(A[1,1]*A[2,3]-A[2,1]*A[1,3])-
A[0,1]*(A[1,0]*A[2,3]-A[2,0]*A[1,3])+
A[0,3]*(A[1,0]*A[2,1]-A[2,0]*A[1,1]));
Result[3,1]:=+det*(A[0,0]*(A[2,1]*A[3,2]-A[3,1]*A[2,2])-
A[0,1]*(A[2,0]*A[3,2]-A[3,0]*A[2,2])+
A[0,2]*(A[2,0]*A[3,1]-A[3,0]*A[2,1]));
Result[3,2]:=-det*(A[0,0]*(A[1,1]*A[3,2]-A[3,1]*A[1,2])-
A[0,1]*(A[1,0]*A[3,2]-A[3,0]*A[1,2])+
A[0,2]*(A[1,0]*A[3,1]-A[3,0]*A[1,1]));
Result[3,3]:=+det*(A[0,0]*(A[1,1]*A[2,2]-A[2,1]*A[1,2])-
A[0,1]*(A[1,0]*A[2,2]-A[2,0]*A[1,2])+
A[0,2]*(A[1,0]*A[2,1]-A[2,0]*A[1,1]));
end;
end;
function Invert(A:TMatrixNM):TMatrixNM;Overload;
var
Sigma:TMatrixNM;
B,U,V:TMatrixNM;
q:TVectorN;
i,N,M:Integer;
begin
SVD(q,u,v,A,True,True);
N:=Length(q);
M:=Length(V);
For i:=0 to N-1 do
begin
q[i]:=1/q[i];
end;
SetLength(Sigma,Length(v),Length(u));
for i:=0 to Length(q)-1 do
Sigma[i,i]:=q[i];
B:=MatrixMulMatrix(Sigma,Transpose(U));
Result:=MatrixMulMatrix(V,B);
end;
-
> Result:=ZeroMatrix44;
Кстати для не до определённых матриц можно так же псевдо обращение применять.
-
Т.е. достаточно 3х3 и поменять знаки у смещения?
Я так сделал. Вроде работает, но медленно:
function InverseMatrix3x3(M: TMatrix4): TMatrix4;
var
detA: float; //Определитель матрицы
NM: TMatrix4; //Матрица миноров
c, r: integer;
begin
//Для отладки
{M[0,0] := 2; M[0,1] := 1; M[0,2] := 1;
M[1,0] := 3; M[1,1] := 2; M[1,2] := 1;
M[2,0] := 2; M[2,1] := 1; M[2,2] := 2;}
//Шаблон матрицы 3×3
// a b c
//[0,0][0,1][0,2] //0
//[1,0][1,1][1,2] //1
//[2,0][2,1][2,2] //2
//1. Находим определитель матрицы
//Формула
//detA = a1*b2*c3 + a3*b1*c2 + a2*b3*c1 - a3*b2*c1 - a1*b3*c2 - a2*b1*c3
//Коррекция индексов массива
//detA = a0*b1*c2 + a2*b0*c1 + a1*b2*c0 - a2*b1*c0 - a0*b2*c1 - a1*b0*c2
detA := M[0,0] * M[1,1] * M[2,2] +
M[0,2] * M[1,0] * M[2,1] +
M[0,1] * M[1,2] * M[2,0] -
M[0,2] * M[1,1] * M[2,0] -
M[0,0] * M[1,2] * M[2,1] -
M[0,1] * M[1,0] * M[2,2];
//Обратная матрица существует при (detA ≠ 0)
if (detA <> 0) then
begin
//2. Находим матрицу миноров
//3. и учитываем матрицу алгебраических дополнений
//[+, -, +]
//[-, +, -]
//[+, -, +]
NM[0,0] := GetMinor3(M, 0, 0);
NM[0,1] := -GetMinor3(M, 0, 1);
NM[0,2] := GetMinor3(M, 0, 2);
NM[1,0] := -GetMinor3(M, 1, 0);
NM[1,1] := GetMinor3(M, 1, 1);
NM[1,2] := -GetMinor3(M, 1, 2);
NM[2,0] := GetMinor3(M, 2, 0);
NM[2,1] := -GetMinor3(M, 2, 1);
NM[2,2] := GetMinor3(M, 2, 2);
//4. Находим транспонированную матрицу алгебраических дополнений
NM := TransponseMatrix3(NM);
//Общий множитель для матрицы
detA := (1 / detA);
//5. Умножаем транспонированную матрицу на определитель
//по формуле 1 / detA * mT
for c := 0 to 2 do //Столбцы
for r := 0 to 2 do //Строки
Result[c, r] := NM[c, r] * detA;
//Проверка для отладки
//NM := MulMatrix(M, Result);
end
else Result := M;
end;
function GetMinor3(var M: TMatrix4; Row, Col: integer): float;
var
r, c: integer;
N: TMatrix2;
NC, NR: integer;
begin
NC := 0; //Строки
NR := 0; //Столбцы
//Столбцы
for c := 0 to 2 do
begin
//Исключаем столбец матрицы
if (c <> Col) then
begin
//Строки
for r := 0 to 2 do
begin
//Исключаем строку матрицы
if (r <> Row) then
begin
N[NR, NC] := M[r, c];
Inc(NC);
end;
end;
Inc(NR);
NC := 0;
end;//if
end;//for j
// a b
//[0,0][0,1] //0
//[1,0][1,1] //1
//Формула: detA = a0*b1 - a1*b0
Result := N[0,0] * N[1,1] - N[1,0] * N[0, 1];
end;
function TransponseMatrix3(M: TMatrix4): TMatrix4;
begin
// a b c
//[0,0][0,1][0,2] //0
//[1,0][1,1][1,2] //1
//[2,0][2,1][2,2] //2
Result[0,0] := M[0,0];
Result[1,0] := M[0,1];
Result[2,0] := M[0,2];
Result[0,1] := M[1,0];
Result[1,1] := M[1,1];
Result[2,1] := M[1,2];
Result[0,2] := M[2,0];
Result[1,2] := M[2,1];
Result[2,2] := M[2,2];
end;
-
В играх матрица по методу гаусса вроде считается. По крайней мере в Quake это так.
-
> Я так сделал. Вроде работает, но медленно:
Раскрутите циклы и заинлайте функции. Хотя детерминант лучше неинлайнить.
function Determinant(const A:TMatrix33):Real;
begin
Result:=(A[0,0]*A[1,1]*A[2,2]+A[0,1]*A[1,2]*A[2,0]+A[0,2]*A[1,0]*A[2,1])-
(A[0,2]*A[1,1]*A[2,0]+A[0,1]*A[1,0]*A[2,2]+A[0,0]*A[1,2]*A[2,1]);
end;
function Invert(A:TMatrix33):TMatrix33;Overload;
var det:real;
begin
det:=Determinant(A);
if det=0 then
begin
Result:=ZeroMatrix33;
end else
begin
det:=1/det;
Result[0,0]:=+det*(A[1,1]*A[2,2]-A[1,2]*A[2,1]);
Result[0,1]:=-det*(A[0,1]*A[2,2]-A[0,2]*A[2,1]);
Result[0,2]:=+det*(A[0,1]*A[1,2]-A[0,2]*A[1,1]);
Result[1,0]:=-det*(A[1,0]*A[2,2]-A[1,2]*A[2,0]);
Result[1,1]:=+det*(A[0,0]*A[2,2]-A[0,2]*A[2,0]);
Result[1,2]:=-det*(A[0,0]*A[1,2]-A[0,2]*A[1,0]);
Result[2,0]:=+det*(A[1,0]*A[2,1]-A[1,1]*A[2,0]);
Result[2,1]:=-det*(A[0,0]*A[2,1]-A[0,1]*A[2,0]);
Result[2,2]:=+det*(A[0,0]*A[1,1]-A[0,1]*A[1,0]);
end;
end;
> В играх матрица по методу гаусса вроде считается. По крайней
> мере в Quake это так.
Quake в плане матриц очень слаб.
> Т.е. достаточно 3х3 и поменять знаки у смещения?
Не совсем. Вечером доберусь до нужного движка напишу. Я просто это случай непроверял. У меня задачи другие, где матрица 4х4 не обязана быть афинной.
-
|
Есть новые 
Нет новых [134427 +38][b:0][p:0]