我通过反投影(Z来自深度缓冲区)获取世界坐标中的鼠标位置。
问题在于这涉及到一个矩阵求逆。当同时使用大和小数值时(例如从原点平移并缩放以查看更多世界),计算变得不稳定。
为了查看这个逆矩阵的精确度,我查看行列式。理想情况下它永远不会为零,因为变换矩阵的本质如此。我知道'det'的小值本身无意义,它可能是由于矩阵中的小值所致。但它也可能是数字出错的迹象。
我还知道可以通过反转每个变换并相乘来计算逆矩阵。这是否提供更高的精度?
我如何判断矩阵是否退化,遭受数值问题?
我如何判断矩阵是否退化,遭受数值问题?
首先,可以查看理解4x4齐次变换矩阵
。Improving accuracy for cumulative matrices (Normalization)
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 orthonormal matrices so no skew or projections ... Orthogonal matrices must be scaled to orthonormal then inverted and then scaled back to make this usable.
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 ...
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
}
//---------------------------------------------------------------------------
The vector operations looks like this:
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];
}
Improving accuracy for non cumulative matrices
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).
Cheap alternative is using double
precision functions like
glTranslated
glRotated
glScaled
...
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 float
s.
Sometimes relative reference frame helps (so keep operations on similar magnitude values) for example see:
ray and ellipsoid intersection accuracy improvement
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.
Pseudo inverse matrix
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 orthonormal 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;
}
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);
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 orthonormal homogenuous 4x4 matrices !!!*
Detection of wrong inverse
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
should be close to 0.0
C
should be close to +1.0
经过一些实验,我发现(谈到变换,不是任意矩阵),矩阵的对角线(即缩放因子)(在反转之前的m
)是行列式值的主要负责者。
因此,我将乘积p = m [0]· m [5]· m [10]· m [15]
(如果它们都!= 0)与行列式进行比较。如果它们相似0.1<p/det<10
,我可以在某种程度上“信任”逆矩阵。否则,我会遇到数值问题,建议更改渲染策略。
near = distance(camera, centerOfWorld) - radusOfWorld
和far = distance(camera, centerOfWorld) + radusOfWorld
。当相机在包围盒内部时,near=nearMin
(例如=1个单位,以查看细节),而far=2*radiusOfWorld
。在这种情况下,我不会关注Z-fighting。 - Ripi2