Всем привет. При изучении 3D в программировании мы пользуемся математикой - линейная алгебра, матрицы, векторы, кватернионы и т.д. В какой-то момент становится интересно как устроено пространство 3D, какие принципы и основы заложены в фундамент отображения моделей. Так же, просто на отображении 3D не возможно остановиться — хочется добавить свет, тень, как минимум. Для расчета света нам надо отправить в шейдер матрицу модели рисуемого объекта — текущего, нормали плоскостей (например на триангулированный объект на каждый треугольник по нормали выходит).
Так же по расчету света необходимо добавить эту же текущую модель, в виде обратной матрицы и транспонированную.
В этой статье покажу как я решил такую задачу - комплексное решение приводящее к обратной матрице (размер 3х3).
Рабочее пространство
Все расчеты буду вести в среде Линукс, 14.2 gcc.
Давайте посмотрим, как можно это решить
Так же подготовим ориентиры которыми будем пользоваться
Так же будем пользоваться транспонированием результирующей матрицы.
Когда мы посчитаем матрицу Minor и применим знаки матрицы Cofactor и применим матрицу Транспонирования - это можно назвать матрицей Сопряжения(Ajoint или она же где-то может называться Ajugate).
Примерно перед применением Транспонирования мы можем проверить определитель (Determinant, Det), определитель нужно проверить на не ноль далее ниже мы рассмотрим случаи 0 и не 0.
Определитель является важной составляющей при расчете обратной матрицы, так как при нахождении обратной матрицы на определитель надо поделить.
Перейдём к реализации
Скрытый текст
#pragma GCC optimize ("Ofast,unroll-loops,-ffast-math")
#pragma GCC target ("sse,sse4.2")
#include <xmmintrin.h>
#include <smmintrin.h>
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <math.h>
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
//обозначим необходимый минимум по структурам
typedef float f32;
struct vec3_t
{
f32 x;
f32 y;
f32 z;
};
typedef struct vec3_t vec3;
struct mat3_t
{
f32 v[9];
};
typedef struct mat3_t mat3;
//функция печати матрицы
void printMatrix(const mat3 m)
{
printf("%f %f %f\n%f %f %f\n%f %f %f\n\n",
m.v[0],m.v[1],m.v[2],
m.v[3],m.v[4],m.v[5],
m.v[6],m.v[7],m.v[8]
);
}
//установка матрицы через 3 вектора
const mat3 setMatrix(const vec3 v,const vec3 v1,const vec3 v2)
{
return (mat3)
{
v.x,v.y,v.z,
v1.x,v1.y,v1.z,
v2.x,v2.y,v2.z,
};
}
//матрица транспонирования
const mat3 Transposedm3(const mat3 m)
{
return (mat3)
{
m.v[0], m.v[3],m.v[6],
m.v[1], m.v[4],m.v[7],
m.v[2], m.v[5],m.v[8]
};
}
//матрица минора
mat3 Minorm3(mat3 m)
{
return (mat3)
{
(m.v[4]*m.v[8]-m.v[5]*m.v[7]),
(m.v[3]*m.v[8]-m.v[5]*m.v[6]),
(m.v[3]*m.v[7]-m.v[4]*m.v[6]),
(m.v[1]*m.v[8]-m.v[2]*m.v[7]),
(m.v[0]*m.v[8]-m.v[2]*m.v[6]),
(m.v[0]*m.v[7]-m.v[1]*m.v[6]),
(m.v[1]*m.v[5]-m.v[2]*m.v[4]),
(m.v[0]*m.v[5]-m.v[2]*m.v[3]),
(m.v[0]*m.v[4]-m.v[1]*m.v[3]),
};
}
//применение знаковой матрицы
mat3 Cofactorm3(mat3 m)
{
return (mat3)
{
m.v[0],-m.v[1],m.v[2],
-m.v[3], m.v[4],-m.v[5],
m.v[6],-m.v[7],m.v[8]
};
}
//результирующая функция Инверсной матрицы
const mat3 Inversem3(mat3 m)
{
mat3 as;
as=Minorm3(m);
as=Cofactorm3(as);
float DetA=(m.v[0]*as.v[0])+(m.v[1]*as.v[1])+(m.v[2]*as.v[2]);
printf("Det %f\n\n",DetA);
//проверка !0
DetA=1/DetA;
as=Transposedm3(as);
return (mat3)
{
as.v[0]*DetA,as.v[1]*DetA,as.v[2]*DetA,
as.v[3]*DetA,as.v[4]*DetA,as.v[5]*DetA,
as.v[6]*DetA,as.v[7]*DetA,as.v[8]*DetA
};
}
int main()
{
//printf("Chapter: Find Inverse Matrix!\nExample Matrix for start calculate:\n\n");
mat3 matA;
matA=setMatrix(
(vec3){1,2,3},
(vec3){4,5,6},
(vec3){7,8,9}
);
printMatrix(matA);
mat3 as;
as=Inversem3(matA);
printf("Inverse Result: \n");
printMatrix(as);
printf("\n\n");
return 0;
}
скомпилируем
gcc -Ofast -ffast-math -o test main.c -lmрезультат
теперь мы умеем находить Обратную матрицу.
Ресурсы
https://www.mathsisfun.com/algebra/matrix-calculator.html
https://www.geeksforgeeks.org/inverse-of-3x3-matrix/
Комментарии (39)
Gay_Lussak
09.12.2024 20:00А вы проверяли алгоритм на устойчивость? Проверьте на матрице Гильберта.
https://ru.wikipedia.org/wiki/Матрица_Гильберта
Jijiki Автор
09.12.2024 20:00спасибо, проверил с сайтом предварительно приведя дроби в десятичные числа,
Скрытый текст
но я перевел дроби в десятичные, у меня и на сайте поидее одинаково
(на обзоре матрица 3х3 простите)
правда сейчас присмотрелся на сайте точность сайта осталась у меня немного округления появились
Jijiki Автор
09.12.2024 20:00(ниже касаемо не нормалей)
библиотека написана вся в таком стиле, пользуюсь ею как глм в преобразованиях неоднократно запускаю визуал, пока только 1 мини проблема, когда двигаю мышку чтобы камера двигалась по вертикали - взор(или угол взора - чтото типо следование направления вектора взгляда) чуть скругляется в центре, это я пока 100 процентов не знаю, в остальном всё как обычно, тоесть в коде матрицы и кватернионы, шейдер считает только три завершающих преобразования
OldFisher
09.12.2024 20:00Если посчитать обратную матрицу аналитически, то правильный ответ получается такой:
9 -36 30 -36 192 -180 30 -180 180Jijiki Автор
09.12.2024 20:00покажите что вы посчитали ?
OldFisher
09.12.2024 20:00Посчитал в Maxima (здесь можно попробовать онлайн http://www.dma.ufv.br/maxima/) такое:
invert(matrix([1, 1/2, 1/3],[1/2, 1/3, 1/4],[1/3, 1/4, 1/5]));Jijiki Автор
09.12.2024 20:00сейчас найду еще какойнибудь ресурс помимо того какой я указал и вы
wataru
09.12.2024 20:00Зачем искать какие-то ресурсы? Если хотите проверить, можно просто перемножить эти 2 матрицы хотя бы вручную.
Jijiki Автор
09.12.2024 20:00понимаю что там перемножить можно чтобы получилась единичная матрица, но ориентиром колебратором выбран результат там где можно посчитать, я таким макаром разобрался почти со всей математикой, метод такой плохой, но был первый неудачный опыт в калибровке отстраивании проверки всей библиотеки и мне тогда было плохо, потомучто там ошибку практически нереально было найти по всей библиотеке, а так я вижу результат пытаюсь понять что не так, тоесть для примера сейчас вся базовая математика работает в визуализаторе(рендере глм при этом намерено выведен из определений) тоесть сейчас проще понять где ошибка. т.к. вся база уже рисует как должно (сейчас первый этап я сравнил результат, убедился что тот калькулятор решает примеры в сторонних ресурсах матриц, разобрался как посчитать, сравнился с калькулятором, второй этап отправить на шейдер матрицу 3х3 и проверить свет)
Jijiki Автор
09.12.2024 20:00вобщем там с 1 ошибкой(тут ошибка это понимание самого принципа по работе с поверхностями(тоесть матрица модель и проче и прочее)). но свет появился (тоесть у новичка даже на глм будет эта ошибка. вобщем я считаю счет работает)(для конкретики для нормализации гладкой поверхности нужно применить Translate затем Scale это минимальный базис для ровной гладкой поверхности чтобы по ней считалась нормаль в шейдере)
Скрытый текст
mat3 t{ modelMatrix.v[0],modelMatrix.v[1],modelMatrix.v[2], modelMatrix.v[4],modelMatrix.v[5],modelMatrix.v[6], modelMatrix.v[8],modelMatrix.v[9],modelMatrix.v[10] }; = Transposedm3(Inversem4(t));
Jijiki Автор
09.12.2024 20:00посмотрел на другом ресурсе да действительно, однако заметим, мой счет отличается от вашего и метод, если вы пользуетесь этим счетом и есть описывающие ресурсы моего метода, расскажите поподробнее что мы считаем с этим калькулятором? и что описывает ресурс по методу расчета тот который я так же указал в ресурсах
пример вывода, все примеры и ход решения в обзоре стороннего ресурса так же были проверены с этим калькулятором
Скрытый текст
OldFisher
09.12.2024 20:00Это не калькулятор, а система компьютерной алгебры. Соответственно, матрицу она обратила не численно, а аналитически, не переходя от дробей к плавающей точки. То есть, решение точное абсолютно и с ним можно сверять правильность своих численных методов. Попутно замечу, что 0.33000 - не очень-то качественное приближение к 1/3, а как раз в случае гильбертовой матрицы это весьма важно.
Jijiki Автор
09.12.2024 20:00Скрытый текст
выводится такое потомучто я посчитал на калькуляторе и подставил значения
Jijiki Автор
09.12.2024 20:00насчет точности наверно, но почему тогда все примеры сторонего ресурса считаются на калькуляторе который я указал в ресурсах, а так же наши расчеты, не в дробях а в десятичных идентичны (кроме округления)?
Gay_Lussak
09.12.2024 20:00Ну у вас очевидно неправильно, вы можете взять сторонний ресурс и убедиться. Собственно я вам и написал, чтобы подтолкнуть на мысль, что ваша реализация через определители не подходит для все матриц.
Jijiki Автор
09.12.2024 20:00нет не очевидно простите, мы используем матрицу транспонирования, подставление алгебраических груп из вводной матрицы и пользуемся матрицей знаков, так же пользуемся определителем, тот метод где получаются другие числа - в его основе лежат совсем другие принципы, например ход решения совсем иной
Jijiki Автор
09.12.2024 20:00значит у задачи есть не 1 ход решения расскажите пожалуйста если вы в этом разбираетесь, потомучто есть калькулятор и есть другие провереннные примеры и сходства есть, и есть метод решения который вы сейчас подразумеваете, и наши результаты не равны в том методе о котором вы говорите и который описал я, я признаться не знал о таком решении не смогу прокоментировать
Jijiki Автор
09.12.2024 20:00простите пожалуйста еще 1 пример
Скрытый текст
Скрытый текст
mat3x3((55.551342, -277.756134, 255.535492), (-277.756134, 1445.919067, -1349.101196), (255.535492, -1349.101196, 1269.742676))Gay_Lussak
09.12.2024 20:00Вам выше уже ответили, вы вместо 1/3 подставляете 0.33. Вы очень грубо округлили. Даже если у вас все будет хорошо после исправления, советую прогнать алгоритм на других плохо-обусловленных матрицах. Если в алгоритме есть деление - это уже повод внимательно присматриваться к нему.
Jijiki Автор
09.12.2024 20:00спасибо тут деление происходит только 1 раз, все проверил и замерил время
1000 000
#include <time.h> int main() { //printf("Chapter: Find Inverse Matrix!\nExample Matrix for start calculate:\n\n"); float GlobalTime=0; srand(time(NULL)); for(int i=0;i<1000000;i++){ struct timespec Ta,Tb; clock_gettime(CLOCK_MONOTONIC_RAW,&Ta); mat3 matA; matA=setMatrix( (vec3){rand()%100,rand()%100,rand()%100}, (vec3){rand()%100,rand()%100,rand()%100}, (vec3){rand()%100,rand()%100,rand()%100} ); //printMatrix(matA); //printf("MinorMatrix then Cofacror:\n\n"); //MinorMatrix mat3 as; as=Inversem3(matA); //printf("Inverse Result: \n"); //printMatrix(as); clock_gettime(CLOCK_MONOTONIC_RAW,&Tb); float localTime=(float)(Tb.tv_nsec-(float)(Ta.tv_nsec)); //printf("%d %f\n",i,localTime); GlobalTime+=localTime; } printf("%f\n",GlobalTime); printf("\n\n"); return 0; }154769280.000000 в наносекундах как я понял (расчет грубый не смотрятся милисекунды) (грубо говоря за 0.2 секунды считается милион обратных матриц в худшем случае 0.4)
тот случай подправил, там всё корректно считается, надо только приводить типы если указываются дроби!!!
Gay_Lussak
09.12.2024 20:00Вы некорректно тестируете. Матрицы которые вы генерируете будут все хорошо обусловлены. Вам необходимо сделать так, чтобы на главной диагонали были одновременно большие и малые числа.
Jijiki Автор
09.12.2024 20:00попробую крайний раз, замер почти такой же при рандоме
matA=setMatrix( (vec3){rand()%10,rand()%100,rand()%100}, (vec3){rand()%100,rand()%10+30,rand()%100}, (vec3){rand()%100,rand()%100,rand()%10+50} );из 20 запусков по 1 милиону 4 отказа тоесть минусовое значение замера времени
Jijiki Автор
09.12.2024 20:00прошу прощения експрессно разобрался в вопросе(посмотрел в WalframAlpha - вы правы)
привожу 2 примера теперь на сайте изза не возможности указания приведения типов отстает результат, но по десятичным он показывает верно,
ниже 2 примера
glm mat3x3( (8.999984, -35.999905, 29.999905), (-35.999905, 191.999496, -179.999527), (29.999905, -179.999527, 179.999557))Скрытый текст
IisNuINu
09.12.2024 20:00видно что самостоятельная работа. Обычно для трёхмерных преобразований используют 4х мерные матрицы, чтобы учесть всякие масштабирования объектов.
OldFisher
09.12.2024 20:00Масштабирование как раз умещается в размерности 3, а вот всякие параллельные переносы и к примеру проецирующее преобразование таки требуют размерности 4 и применения так называемых "однородных координат".
IisNuINu
09.12.2024 20:00И вот ещё что, вы простите, поленились, в структуре mat3_t вы задали коэффициэнты матрицы через вектор. это, честно говоря какой то нонсенс. структуры для того и создаются что бы именовать свои поля, компилятор потом всё это корректно обработает. Поэтому если будете переделывать проименуйте отдельно каждый коэффицэнт как нибудь так mIJ, где I номер строки J номер столбца, или как нибудь подобно. Вам потом самому будет проще разбираться с алгоритмом.
Jijiki Автор
09.12.2024 20:00спасибо, о чем вы? не надо этого создавать компилятор ничего не упростит в таком случае посмотрите ассемблер этого кода, касаемо идеи flat пошустрее работает, не только в матрица но и в кубоматематике - кубоматематика продолжение матрицы, где единица центроид вокселя
Jijiki Автор
09.12.2024 20:00вот еще рабочий пример
Скрытый текст
mat4 AngleAxisfv3m4(mat4 m1,float deg,vec3 axis) { float rad = deg * DPR1; if(MagnitudeSqv3(axis)!=EPSILON1) { axis = Normalizev3(axis); } Quaternion res; res.v[0]= axis.x * sinf(rad * 0.5f); res.v[1]= axis.y * sinf(rad * 0.5f); res.v[2]= axis.z * sinf(rad * 0.5f); res.v[3]= cosf(rad * 0.5f); f32 r00=2 * (res.v[0] * res.v[0] + res.v[1] * res.v[1])-1; f32 r01=2 * (res.v[1] * res.v[2] - res.v[0] * res.v[3]); f32 r02=2 * (res.v[1] * res.v[3] + res.v[0] * res.v[2]); f32 r10=2 * (res.v[1] * res.v[2] + res.v[0] * res.v[3]); f32 r11=2 * (res.v[0] * res.v[0] + res.v[2] * res.v[2])-1; f32 r12=2 * (res.v[2] * res.v[3] - res.v[0] * res.v[1]); f32 r20=2 * (res.v[1] * res.v[3] - res.v[0] * res.v[2]); f32 r21=2 * (res.v[2] * res.v[3] + res.v[0] * res.v[1]); f32 r22=2 * (res.v[0] * res.v[0] + res.v[3] * res.v[3])-1; mat4 t=(mat4){ r00,r01,r02,0, r10,r11,r12,0, r20,r21,r22,0, 0,0,0,1 }; return Mulm4(m1,t); }возможно нету только указываюещего ориентирования, но работает даже так
Jijiki Автор
09.12.2024 20:00Скрытый текст
float Mat2ValueAt(mat2 input, int row, int column) { return input.v[column * 2 + row]; } float Mat3ValueAt(mat3 input, int row, int column) { return input.v[column * 3 + row]; } float Mat4ValueAt(mat4 input, int row, int column) { return input.v[column * 4 + row]; }вот еще пример
OldFisher
09.12.2024 20:00Отметим, что если не использовать масштабирование (особенно неравномерное) для размещения объектов в сцене, остаются только повороты и смещения. В описывающей такое размещение матрице (она и называется матрица модели) размером 4*4 смещения будут отдельно в 4-м столбце, а повороты займут подматрицу 3*3 (столбцы и строки с 1 по 3). Если рассмотрим её как отдельную матрицу 3*3, то для получения обратной её достаточно транспонировать. Ещё одно занятное свойство: длина векторов, составляющих её строки и столбцы равна 1.
В итоге, обращение полной матрицы с такими ограничениями сводится к одному транспонированию и одному умножению вектора на матрицу для получения обратных смещений.
Если было задействовано равномерное масштабирование, длина векторов-строк и векторов-столбцов будет равна коэффициенту масштабирования, и соответственно для обращения подматрицу придётся не только транспонировать, но и делить на квадрат коэффициента.
Такое использование частных случаев получается дешевле, и может быть иногда точнее, чем обращать в общем виде.
Jijiki Автор
09.12.2024 20:00если вы про Scale и замечание общее про сцену, то я это использую, под обьект позиция - поворот- скейл, далее их перемножение матричное, а шейдер доводит проекцию модель вид, этот расчет, который представлен производится для матрицы нормалей применяется в расчете света, я посмотрел пример как делают, и так как воспроизвёл всю библиотеку не было такого метода чтобы найти принципиально математически обратную матрицу
, единственное что не понял покачто что происходит с матрицей в том же глм, передаем в матрицу размера 3 матрицу 4 в шейдере из матрицы получается вектор, и свет работает, вот интересно 4 поле просто отсекается или там какойто расчет
Jijiki Автор
09.12.2024 20:00тот пример что я смотрел, конкретно выглядит так
mat4 ExampleObject;
на входе у шейдера mat3 ; glm::transpose(glm::inverse(ExampleObject));
а в шейдере он домножая на вектор получает вектор -> вектор пошел в свет
Jijiki Автор
09.12.2024 20:00добавлю еще всё так просто у меня OpenGL весь математический визуал подведён под него если надо могу кинуть визуал, везде написано про 4 столбец но у меня почемуто 4 строка
Скрытый текст
mat4 setOrthoFrustumm4(const float l,const float r,const float b,const float t,const float n,const float f) { return (mat4){ 2.0f / (r - l),0,0, 0, 0, 2.0f / (b - t), 0, 0, 0, 0, -2.0f / (f - n),0, -(r + l) / (r - l),-(t + b) / (b - t),-(f + n) / (f - n),1 }; }вот яркий пример, смотрел обзор как считают по математике со столбцом я так и не понял а тут всё заработало как надо, возможно вся суть держится на перемножающей матрице, ну и признаюсь у меня всё в векторах хранится вроде, хотя помню когда вооброжал гдето поидее проскальзывает переворот-переворот хранения данных со строк на столбцы, щас не вспомню где именно
Jijiki Автор
09.12.2024 20:00проверил, что вы написали и с обратной матрицей и с просто транспонированием 1 и тоже (в свете)
mat3 t{ modelMatrix.v[0],modelMatrix.v[1],modelMatrix.v[2], modelMatrix.v[4],modelMatrix.v[5],modelMatrix.v[6], modelMatrix.v[8],modelMatrix.v[9],modelMatrix.v[10] }; Transposedm3(t); mat3 t{ modelMatrix.v[0],modelMatrix.v[1],modelMatrix.v[2], modelMatrix.v[4],modelMatrix.v[5],modelMatrix.v[6], modelMatrix.v[8],modelMatrix.v[9],modelMatrix.v[10] }; = Transposedm3(Inversem4(t));просто одно и тоже, спасибо вам, ну так хоть разобрался вроде досконально что и как
понял что тут много хитростей по оптимизации(что как и когда умножать - преобразовывать)
Jetstorm
Хоть кто-то в рукопашную делает, а не использует glm или что-то подобное
Jijiki Автор
советую, если заинтересованы, интересно достаточно, я почти всю библиотеку реализовал
сверяюсь по сайту который прикреплен в ресурсах, по нему почти всё что нужно проверил, изучил