Матрица обратной точности

Question

У меня большой мир, около 5 000 000 х 1 000 000 единиц. Камера может находиться рядом с каким-либо объектом или достаточно далеко, чтобы видеть весь мир.
Я получаю положение мыши в мировых координатах, не проектируя (Z приходит из буфера глубины). Проблема в том, что она включает в себя обратную матрицу. При одновременном использовании больших и малых чисел (например, перевод от исходного положения и масштабирование, чтобы увидеть больше мира), вычисления становятся нестабильными.

Пытаясь увидеть точность этой обратной матрицы, я смотрю на определитель. В идеале оно никогда не будет нулевым из-за природы матриц преобразования. Я знаю, что "det" небольшое значение само по себе ничего не значит, это может быть связано с небольшими значениями в матрице. Но это также может быть признаком того, что цифры становятся неправильными.

Я также знаю, что могу вычислить обратное, инвертируя каждое преобразование и умножая их. Это обеспечивает большую точность?

Как я могу сказать, что моя матрица дегенерирует, страдают от численных проблем?

Answer 1

Для начала см. Понимание матриц однородного преобразования 4x4

1. Повышение точности для кумулятивных матриц (нормализация)

   Чтобы избежать вырождения матрицы преобразования, выберите одну ось в качестве главной. Я обычно выбирал Z как это обычно вид или направление вперед в моих приложениях. Затем используйте перекрестное произведение для пересчета / нормализации остальных осей (которые должны быть перпендикулярны друг другу и, если не используется масштаб, то также размер единицы). Это может быть сделано только для ортогональных / ортонормированных матриц, поэтому нет перекоса или проекций...

   Вам не нужно делать это после каждой операции, просто сделайте счетчик операций, выполненных на каждой матрице, и, если какой-то порог пройден, нормализуйте его и сбросьте счетчик.

   Чтобы обнаружить вырождение таких матриц, вы можете проверить ортогональность по точечному произведению между любыми двумя осями (должно быть ноль или очень близко к нему). Для ортонормированных матриц вы можете проверить также размер блока векторов направления оси...

   Вот как выглядит нормализация моей матрицы преобразования (для ортонормированных матриц) в C++:

   double reper::rep[16]; // this is my transform matrix stored as member in `reper` class
   //---------------------------------------------------------------------------
   void reper::orto(int test) // test is for overiding operation counter
           {
           double   x[3],y[3],z[3]; // space for axis direction vectors
           if ((cnt>=_reper_max_cnt)||(test)) // if operations count reached or overide
                   {
                   axisx_get(x);      // obtain axis direction vectors from matrix
                   axisy_get(y);
                   axisz_get(z);
                   vector_one(z,z);   // Z = Z / |z|
                   vector_mul(x,y,z); // X = Y x Z  ... perpendicular to y,z
                   vector_one(x,x);   // X = X / |X|
                   vector_mul(y,z,x); // Y = Z x X  ... perpendicular to z,x
                   vector_one(y,y);   // Y = Y / |Y|
                   axisx_set(x);      // copy new axis vectors into matrix
                   axisy_set(y);
                   axisz_set(z);
                   cnt=0;             // reset operation counter
                   }
           }
   
   //---------------------------------------------------------------------------
   void reper::axisx_get(double *p)
           {
           p[0]=rep[0];
           p[1]=rep[1];
           p[2]=rep[2];
           }
   //---------------------------------------------------------------------------
   void reper::axisx_set(double *p)
           {
           rep[0]=p[0];
           rep[1]=p[1];
           rep[2]=p[2];
           cnt=_reper_max_cnt; // pend normalize in next operation that needs it
           }
   //---------------------------------------------------------------------------
   void reper::axisy_get(double *p)
           {
           p[0]=rep[4];
           p[1]=rep[5];
           p[2]=rep[6];
           }
   //---------------------------------------------------------------------------
   void reper::axisy_set(double *p)
           {
           rep[4]=p[0];
           rep[5]=p[1];
           rep[6]=p[2];
           cnt=_reper_max_cnt; // pend normalize in next operation that needs it
           }
   //---------------------------------------------------------------------------
   void reper::axisz_get(double *p)
           {
           p[0]=rep[ 8];
           p[1]=rep[ 9];
           p[2]=rep[10];
           }
   //---------------------------------------------------------------------------
   void reper::axisz_set(double *p)
           {
           rep[ 8]=p[0];
           rep[ 9]=p[1];
           rep[10]=p[2];
           cnt=_reper_max_cnt; // pend normalize in next operation that needs it
           }
   //---------------------------------------------------------------------------
   

   Векторные операции выглядят так:

   void  vector_one(double *c,double *a)
           {
           double l=divide(1.0,sqrt((a[0]*a[0])+(a[1]*a[1])+(a[2]*a[2])));
           c[0]=a[0]*l;
           c[1]=a[1]*l;
           c[2]=a[2]*l;
           }
   void  vector_mul(double *c,double *a,double *b)
           {
           double   q[3];
           q[0]=(a[1]*b[2])-(a[2]*b[1]);
           q[1]=(a[2]*b[0])-(a[0]*b[2]);
           q[2]=(a[0]*b[1])-(a[1]*b[0]);
           for(int i=0;i<3;i++) c[i]=q[i];
           }
   

2. Повышение точности для не кумулятивных матриц

   Ваш единственный выбор - использовать хотя бы double Точность ваших матриц. Безопаснее всего использовать GLM или собственную матричную математику, основанную хотя бы на double тип данных (как мой reper учебный класс).

   Дешевая альтернатива использует double точные функции, такие как

   glTranslated
   glRotated
   glScaled
   ...
   

   что в некоторых случаях помогает, но не безопасно, так как реализация OpenGL может усекать его до float, Также еще нет 64-битных HW- интерполяторов, поэтому все повторяющиеся результаты между этапами конвейера усекаются до floats.

   Иногда помогает относительная система отсчета (поэтому сохраняйте операции с одинаковыми значениями величины), например, смотрите:

   • повышение точности пересечения лучей и эллипсоидов

   Кроме того, в случае, если вы используете собственные математические математические функции, вы должны учитывать также порядок операций, чтобы вы всегда теряли наименьшую возможную точность.

3. Псевдообратная матрица

   В некоторых случаях вы можете избежать вычисления обратной матрицы с помощью детерминантов, схемы Хорнера или метода исключения Гаусса, потому что в некоторых случаях вы можете использовать тот факт, что Transpose ортогональной вращательной матрицы также является ее обратной. Вот как это делается:

   void  matrix_inv(GLfloat *a,GLfloat *b) // a[16] = Inverse(b[16])
           {
           GLfloat x,y,z;
           // transpose of rotation matrix
           a[ 0]=b[ 0];
           a[ 5]=b[ 5];
           a[10]=b[10];
           x=b[1]; a[1]=b[4]; a[4]=x;
           x=b[2]; a[2]=b[8]; a[8]=x;
           x=b[6]; a[6]=b[9]; a[9]=x;
           // copy projection part
           a[ 3]=b[ 3];
           a[ 7]=b[ 7];
           a[11]=b[11];
           a[15]=b[15];
           // convert origin: new_pos = - new_rotation_matrix * old_pos
           x=(a[ 0]*b[12])+(a[ 4]*b[13])+(a[ 8]*b[14]);
           y=(a[ 1]*b[12])+(a[ 5]*b[13])+(a[ 9]*b[14]);
           z=(a[ 2]*b[12])+(a[ 6]*b[13])+(a[10]*b[14]);
           a[12]=-x;
           a[13]=-y;
           a[14]=-z;
           }
   

   Таким образом, вращающаяся часть матрицы транспонируется, проекция остается прежней, а исходная позиция пересчитывается A*inverse(A)=unit_matrix Эта функция написана так, что она может быть использована как на месте, так что вызов

   GLfloat a[16]={values,...}
   matrix_inv(a,a);
   

   привести к действительным результатам тоже. Этот способ вычисления Inverse является более быстрым и безопасным в численном отношении, поскольку он требует гораздо меньше операций (без рекурсии или сокращения без делений). Конечно, это работает только для ортогональных однородных матриц 4х4!!!*

4. Обнаружение неправильного обратного

   Так что, если у вас есть матрица A и его обратное B затем:

   A*B = C = ~unit_matrix
   

   Так что умножьте обе матрицы и проверьте на единицу матрицы...

   • абсолютная сумма всех недиагональных элементов C должно быть близко к 0.0
   • все диагональные элементы C должно быть близко к +1.0

Answer 2

После некоторых экспериментов я вижу, что (говоря о преобразованиях, а не о какой-либо матрице) диагональ (то есть коэффициенты масштабирования) матрицы (m, до инвертирования) является основной ответственной за определяющее значение.

Поэтому я сравниваю продукт p= m[0] · m[5] · m[10] · m[15] (если все они!= 0) с определителем. Если они похожи 0.1 < p/det < 10 Я могу как-то "доверять" в обратной матрице. В противном случае у меня есть числовые проблемы, которые советуют изменить стратегию рендеринга.