opengl - Precisión inversa de matriz

math matrix

Question

Tengo un mundo grande, alrededor de 5,000,000 x 1,000,000 de unidades. La cámara puede estar cerca de algún objeto o lo suficientemente lejos como para ver el mundo entero.
Obtengo la posición del mouse en coordenadas mundiales desproyectando (Z viene del buffer de profundidad). El problema es que involucra una matriz inversa . Cuando se usan números grandes y pequeños (p. Ej., Traducir desde el origen y escalar para ver más mundo) al mismo tiempo, los cálculos se vuelven inestables.

Tratando de ver la precisión de esta matriz inversa, miro el determinante. Idealmente, nunca será cero, debido a la naturaleza de las matrices de transformación. Sé que ser ''det'' un valor pequeño no significa nada por sí solo, puede deberse a valores pequeños en la matriz. Pero también puede ser una señal de que los números se están equivocando.

También sé que puedo calcular el inverso invirtiendo cada transformación y multiplicándolas. ¿Proporciona más precisión?

¿Cómo puedo saber si mi matriz se está degenerando, sufro problemas numéricos?

Answer 1

Después de algunos experimentos, veo que (hablando de transformaciones, no de ninguna matriz) la diagonal (es decir, los factores de escala) de la matriz ( m , antes de invertir) es el principal responsable del valor determinante.

Así que comparo el producto p= m[0] · m[5] · m[10] · m[15] (si todos son! = 0) con el determinante. Si son similares 0.1 < p/det < 10 puedo "confiar" de alguna manera en la matriz inversa. De lo contrario, tengo problemas numéricos que aconsejan cambiar la estrategia de representación.

Answer 2

para empezar, ver Comprender matrices de transformación homogéneas 4x4

1. Mejora de la precisión para matrices acumulativas (Normalización)

   Para evitar la degeneración de la matriz de transformación, seleccione un eje como principal. Por lo general, elijo Z ya que generalmente es la vista hacia adelante o hacia adelante en mis aplicaciones. Luego explote el producto cruzado para recalcular / normalizar el resto de ejes (que deben ser perpendiculares entre sí y, a menos que se use escala, también el tamaño de la unidad). Esto se puede hacer solo para matrices ortogonales / ortonormales, por lo que no hay sesgo ni proyecciones ...

   No necesita hacer esto después de cada operación, solo haga un contador de operaciones realizadas en cada matriz y si se cruza algún umbral, normalícelo y reinicie el contador.

   Para detectar la degeneración de tales matrices, puede probar la ortogonalidad por producto de punto entre dos ejes (debe ser cero o muy cerca). Para matrices ortonormales, puede probar también el tamaño unitario de los vectores de dirección del eje ...

   Así es como se ve mi normalización de matriz de transformación (para matrices ortonormales ) en 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
           }
   //---------------------------------------------------------------------------

   Las operaciones vectoriales se ven así:

   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. Mejora de la precisión para matrices no acumulativas

   Su única opción es utilizar al menos el double precisión de sus matrices. Lo más seguro es usar GLM o su propia matriz matemática basada al menos en double tipo de datos double (como mi clase de reper ).

   Una alternativa barata es usar funciones de double precisión como

   glTranslated
   glRotated
   glScaled
   ...

   lo que en algunos casos ayuda pero no es seguro ya que la implementación de OpenGL puede truncarlo para que float . Además, todavía no hay interpoladores HW de 64 bits, por lo que todos los resultados iterados entre las etapas de la tubería se truncan a float .

   A veces, el marco de referencia relativo ayuda (así que mantenga las operaciones en valores de magnitud similares), por ejemplo, vea:

   • mejora de la precisión de la intersección de rayos y elipsoides

   Además, en caso de que esté utilizando sus propias funciones matemáticas de matriz, debe tener en cuenta también el orden de las operaciones para que siempre pierda la menor cantidad de precisión posible.

3. Matriz pseudo inversa

   En algunos casos, puede evitar el cálculo de la matriz inversa mediante determinantes o el esquema de Horner o el método de eliminación de Gauss porque en algunos casos puede explotar el hecho de que la Transposición de la matriz rotacional ortogonal también es inversa . Así es como se hace:

   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;
           }

   Por lo tanto, la parte rotacional de la matriz se transpone, la proyección se mantiene como estaba y la posición de origen se vuelve a calcular de manera que A*inverse(A)=unit_matrix Esta función se escribe para que se pueda usar como in situ.

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

   conducir a resultados válidos también. Esta forma de calcular Inverse es más rápida y numéricamente más segura, ya que requiere muchas menos operaciones (sin recursiones ni reducciones, sin divisiones ). ¡De grosor, esto funciona solo para matrices ortogonales homogéneas de 4x4! *

4. Detección de inversa incorrecta

   Entonces, si tienes la matriz A y su inversa B entonces:

   A*B = C = ~unit_matrix

   Entonces multiplique ambas matrices y verifique la matriz unitaria ...

   • La suma absoluta de todos los elementos no diagonales de C debe estar cerca de 0.0
   • todos los elementos diagonales de C deben estar cerca de +1.0