Мастера DELPHI, Delphi programming community Рейтинг@Mail.ru Титульная страница Поиск, карта сайта Написать письмо 
| Новости |
Новости сайта
Поиск |
Поиск по лучшим сайтам о Delphi
FAQ |
Огромная база часто задаваемых вопросов и, конечно же, ответы к ним ;)
Статьи |
Подборка статей на самые разные темы. Все о DELPHI
Книги |
Новинки книжного рынка
Новости VCL
Обзор свежих компонент со всего мира, по-русски!
|
| Форумы
Здесь вы можете задать свой вопрос и наверняка получите ответ
| ЧАТ |
Место для общения :)
Орешник |
Коллекция курьезных вопросов из форумов
KOL и MCK |
KOL и MCK - Компактные программы на Delphi
Основная («Начинающим»)/ Базы / WinAPI / Компоненты / Сети / Media / Игры / Corba и COM / KOL / FreePascal / .Net / Прочее / rsdn.org

 
Чтобы не потерять эту дискуссию, сделайте закладку « предыдущая ветвь | форум | следующая ветвь »

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


dmk ©   (17.05.17 17:11

Все привет! Кто-нибудь знает как инвертировать матрицу 4x3?
3x3  сделал, а 4x3 никак. У меня по методу Крамера. Формул для 4x4 ни у кого не завалялось?


Pavia ©   (18.05.17 08:51[1]

Метод Крамера работает только с однозначно определёнными матрицами.
А матрица 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;


Pavia ©   (18.05.17 08:57[2]


> Result:=ZeroMatrix44;

Кстати для не до определённых матриц можно так же псевдо обращение применять.


dmk ©   (18.05.17 21:31[3]

Т.е. достаточно 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;


dmk ©   (18.05.17 21:33[4]

В играх матрица по методу гаусса вроде считается. По крайней мере в Quake это так.


Pavia ©   (19.05.17 08:34[5]


> Я так сделал. Вроде работает, но медленно:

Раскрутите циклы и заинлайте функции. Хотя детерминант лучше неинлайнить.

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 не обязана быть афинной.


Pavia ©   (19.05.17 10:13[6]

https://github.com/OGRECave/ogre/blob/ac15c70633fd690f85b46411d857eed56e9ed2bb/OgreMain/src/OgreMatrix4.cpp#L164
Там не просто знак меняется но и ещё умножается на результирующую подматрицу 3х3


версия для печати

Написать ответ

Ваше имя (регистрация  E-mail 







Разрешается использование тегов форматирования текста:
<b>жирный</b> <i>наклонный</i> <u>подчеркнутый</u>,
а для выделения текста программ, используйте <code> ... </code>
и не забывайте закрывать теги! </b></i></u></code> :)


Наверх

  Рейтинг@Mail.ru     Титульная страница Поиск, карта сайта Написать письмо