矩阵逆精度

Question

我有一个很大的世界，大约有 5,000,000 x 1,000,000 个单位。相机可以靠近某个物体，也可以远到足以看到整个世界。
我通过不投影（Z 来自深度缓冲区）获得世界坐标中的鼠标位置。 问题是它涉及到一个逆矩阵。当同时使用大小数字（例如从原点平移和缩放以查看更多世界）时，计算变得不稳定。

试图查看此逆矩阵的准确性，我查看了行列式。理想情况下，由于变换矩阵的性质，它永远不会为零。我知道'det'一个小值本身并不意味着什么，这可能是由于矩阵中的小值。但这也可能是数字出错的迹象。

我也知道我可以通过反转每个变换并将它们相乘来计算逆。它是否提供更高的准确性？

我如何判断我的矩阵是否正在退化，是否出现数值问题？

Answer 1

对于初学者，请参阅Understanding 4x4 homogenous transform matrices

1. 提高累积矩阵的准确性（归一化）

   为了避免变换矩阵的退化，选择一个轴作为主轴。我通常选择Z，因为它通常是我的应用程序中的视图或前进方向。然后利用 cross product 重新计算/归一化其余的轴（它们应该相互垂直，除非使用比例尺，否则还有单位大小）。这只能用于正交矩阵，因此不会出现偏斜或投影......正交矩阵必须缩放为正交矩阵，然后反转，然后再缩小以使其可用。

   您不需要在每次操作后执行此操作，只需对每个矩阵进行操作计数器，如果超过某个阈值，则对其进行归一化并重置计数器。

   要检测此类矩阵的退化，您可以通过任意两个轴（应该为零或非常接近）之间的点积来测试正交性。对于正交矩阵，您还可以测试轴方向向量的单位大小...

   这是我的变换矩阵归一化在 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。

   有时相对参考系会有所帮助（因此请保持对相似幅度值的操作），例如：

   ray and ellipsoid intersection accuracy improvement

   此外，如果您使用自己的矩阵数学函数，您还必须考虑运算顺序，这样您总是会损失尽可能少的准确度。

3. 伪逆矩阵

   在某些情况下，您可以避免通过行列式或霍纳方案或高斯消元法计算逆矩阵，因为在某些情况下，您可以利用 正交旋转矩阵的转置也是它的逆矩阵这一事实。以下是它的完成方式：

   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 的方法更快且在数值上更安全，因为它需要更少的操作（没有递归或归约没有除法）。粗略 这仅适用于正交同质 4x4 矩阵！！！*

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，我可以以某种方式“信任”逆矩阵。否则，我会遇到建议更改渲染策略的数值问题。