矩阵逆精度 [英] Matrix inverse accuracy
问题描述
我的世界很大，大约有5,000,000 x 1,000,000单位.相机可以靠近某个物体，也可以足够远以看到整个世界.
我通过不投影获得鼠标在世界坐标中的位置(Z来自深度缓冲区).
问题在于它涉及到矩阵逆.当同时使用大小数字时(例如，远离原点平移并缩放以查看更多世界)，计算将变得不稳定.
试图查看这个逆矩阵的准确性，我看了行列式.理想情况下，由于转换矩阵的性质，它永远不会为零.我知道被删除"一个小值本身并不意味着什么，这可能是由于矩阵中的小值所致.但这也可能是数字出错的迹象.
我也知道我可以通过对每个变换求逆并将它们相乘来计算逆.它是否提供更高的准确性?
我怎么知道我的矩阵是否退化，遇到数值问题?
I have a large world, about 5,000,000 x 1,000,000 units. The camera can be near some object or far enough as to see the whole world.
I get the mouse position in world coordinates by unprojecting (Z comes from depth buffer).
The problem is that it involves a matrix inverse. When using big and small numbers (e.g. translating away from origin and scaling to see more world) at the same time, the calculations become unstable.
Trying to see the accuracy of this inverse matrix I look at the determinant. Ideally it never will be zero, due to the nature of transformation matrices. I know that being 'det' a small value means nothing on its own, it can be due to small values in the matrix. But it can also be a sign of numbers becoming wrong.
I also know I can calculate the inverse by inverting each transformation and multiplicaying them. Does it provide more accuracy?
How can I tell if my matrix is getting degenerated, suffer numerical issues?
推荐答案
对于初学者，请参见了解4x4均匀变换矩阵
-
提高累积矩阵的准确性(归一化)
要避免退化，请选择一个轴作为主轴.我通常选择Z，因为它通常是我的应用程序中的视图或前进方向.然后利用叉积重新计算/归一化其余的轴(轴应彼此垂直，除非使用比例尺，否则还应使用单位尺寸).只能对正交/正交矩阵执行此操作，因此不会产生偏斜或投影...
To avoid degeneration of transform matrix select one axis as main. I usually chose Z as it is usually view or forward direction in my apps. Then exploit cross product to recompute/normalize the rest of axises (which should be perpendicular to each other and unless scale is used then also unit size). This can be done only for orthogonal/orthonormal matrices so no skew or projections ...
您不必在每次操作后都要做此操作，只需对每个矩阵上完成的操作做一个计数器即可，如果超过某个阈值，则将其标准化并重置计数器.
You do not need to do this after every operation just make a counter of operations done on each matrix and if some threshold crossed then normalize it and reset counter.
要检测这种矩阵的退化，您可以通过任意两个轴(应该为零或非常接近)之间的点积来测试正交性.对于正交矩阵，您还可以测试轴方向矢量的单位大小...
To detect degeneration of such matrices you can test for orthogonality by dot product between any two axises (should be zero or very near it). For orthonormal matrices you can test also for unit size of axis direction vectors ...
这是在 C ++ 中我的变换矩阵归一化的外观(对于正交矩阵):
Here is how my transform matrix normalization looks like (for orthonormal matrices) in 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];
}
提高非累积矩阵的准确性
您唯一的选择是至少使用矩阵的double精度.最安全的方法是使用 GLM 或至少基于double数据类型(例如我的reper类)的自己的矩阵数学.
Your only choice is use at least double accuracy of your matrices. Safest is to use GLM or your own matrix math based at least on double data type (like my reper class).
便宜的替代品使用的是double精度函数，例如
Cheap alternative is using double precision functions like
glTranslated
glRotated
glScaled
...
在某些情况下有帮助，但不安全，因为 OpenGL 实现可以将其截断为float.此外，还没有64位 HW 插值器，因此，流水线级之间的所有迭代结果都将被截断为float s.
which in some cases helps but is not safe as OpenGL implementation can truncate it to float. Also there are no 64 bit HW interpolators yet so all iterated results between pipeline stages are truncated to floats.
有时，相对参考框架会有所帮助(因此，请对相似的幅度值进行操作)，例如:
Sometimes relative reference frame helps (so keep operations on similar magnitude values) for example see:
此外，如果您使用自己的矩阵数学函数，则还必须考虑运算顺序，这样一来，您始终可能会损失尽可能小的精度.
Also In case you are using own matrix math functions you have to consider also the order of operations so you always lose smallest amount of accuracy possible.
伪逆矩阵
在某些情况下，您可以避免通过行列式或Horner方案或高斯消除方法来计算逆矩阵，因为在某些情况下，您可以利用正交旋转矩阵的转置也是其逆值这一事实.这是完成的方式:
In some cases you can avoid computing of inverse matrix by determinants or Horner scheme or Gauss elimination method because in some cases you can exploit the fact that Transpose of orthogonal rotational matrix is also its inverse. Here is how it is done:
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编写此函数以便可以就地使用它来调用
So rotational part of the matrix is transposed, projection stays as was and origin position is recomputed so A*inverse(A)=unit_matrix This function is written so it can be used as in-place so calling
GLfloat a[16]={values,...}
matrix_inv(a,a);
也导致有效的结果.这种计算Inverse的方法更快，在数值上更安全，因为它减少了很多操作(无需递归或归约无除法).粗糙的仅适用于正交同构4x4矩阵!!! *
lead to valid results too. This way of computing Inverse is quicker and numerically safer as it pends much less operations (no recursion or reductions no divisions). Of coarse this works only for orthogonal homogenuous 4x4 matrices !!!*
错误的逆向检测
因此，如果您得到矩阵A及其逆矩阵B，那么:
So if you got matrix A and its inverse B then:
A*B = C = ~unit_matrix
因此将两个矩阵相乘并检查单位矩阵...
So multiply both matrices and check for unit matrix...
-
C的所有非对角元素的绝对和应接近0.0 -
C的所有对角线元素应接近+1.0
- abs sum of all non diagonal elements of
Cshould be close to0.0 - all diagonal elements of
Cshould be close to+1.0
这篇关于矩阵逆精度的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
