Exactitud inversa de matriz

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 las coordenadas mundiales al no proyectar (Z proviene del buffer de profundidad). El problema es que implica una matriz inversa . Al usar números grandes y pequeños (por ejemplo, traducir lejos del 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, observo 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 en sí mismo, puede ser debido a pequeños valores en la matriz. Pero también puede ser un signo 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, sufre problemas numéricos?

para empezar, véase Understanding 4×4 homogeneus transform matrices

  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 por lo general se ve o reenvía en mis aplicaciones. Luego explote el producto cruzado para volver a calcular / normalizar el rest de ejes (que deben ser perpendiculares entre sí y, a menos que se use una escala, también el tamaño de la unidad). Esto se puede hacer solo para matrices ortogonales / ortonormales, por lo que no hay desviaciones ni proyecciones …

    No es necesario que haga esto después de cada operación para hacer un contador de operaciones en cada matriz y si se cruza un 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 cualquiera (debe ser cero o muy cercano). Para las matrices ortonormales también puede probar el tamaño unitario de los vectores de dirección del eje …

    Así es como se ve la normalización de la 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 usar al menos double precisión double de sus matrices. Lo más seguro es utilizar GLM o su propia matriz matemática basada, al menos, en el tipo de datos double (como mi clase reper ).

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

     glTranslated glRotated glScaled ... 

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

    A veces, el marco de referencia relativo ayuda (por lo tanto, mantenga las operaciones en valores de magnitud similares), por ejemplo, consulte:

    • mejora de la precisión de intersección de rayo y elipsoide

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

  3. Pseudo matriz inversa

    En algunos casos puede evitar la computación de la matriz inversa por 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 de rotación 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; } 

    De modo que la parte rotatoria de la matriz se transpone, la proyección se mantiene como estaba y la posición de origen se vuelve a calcular de modo que A*inverse(A)=unit_matrix Esta función se escribe para que pueda usarse en el lugar de modo que llame

     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 menos operaciones (sin recursión ni reducciones ni divisiones ). ¡De grosor, esto solo funciona para matrices ortogonales 4 × 4 homogenuas! *

  4. Detección de inversa incorrecta

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

     A*B = C = ~unit_matrix 

    Así que multiplique ambas matrices y verifique la matriz de la unidad ...

    • la sum abs de todos los elementos no diagonales de C debería estar cerca de 0.0
    • todos los elementos diagonales de C deberían estar cerca de +1.0

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 la principal responsable del valor determinante.

Entonces comparo el producto p= m[0] · m[5] · m[10] · m[15] (si todos ellos 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 para el renderizado.